demo_fem2.m

Go to the documentation of this file.

function demo_fem2()
%function demo_fem2()
% Demo of +Fem usage for finite elements simulations of linear and
% stationary models: (i) Poisson, (ii) Stokes and nonlinear models:
% (iii) Navier-Stokes.

% IM, 09.09.2016

% set up model:
fprintf('\n');
disp('(a) P_k-FEM on triangular grid for Poisson problem');
disp('==================================================');
params = [];
params.pdeg = 2;
params.numintervals = 2*2^4;
model = poisson_model(params);
model = generic_fem_model_adapter(model);
model_data = gen_model_data(model);

% Simulation:
sim_data = detailed_simulation(model, model_data);

% plot:
disp(['ndofs = ',num2str(sim_data.uh.ndofs)]);
disp('plotting solution ...');
figure, plot_sim_data(model, model_data, sim_data, []);
title('FEM solution of -Laplace u = f for zero-boundary values')

disp('press enter to continue');
pause;
close all

% eoc:
fprintf('\n');
disp('error convergence table:')

i_s = [0,1,2,3,4];
nis = length(i_s);
l2errs = zeros(nis,1);
h1errs = zeros(nis,1);
ndofs = zeros(nis,1);

for pdeg = 1:2
    fprintf('\n');
    disp(['              *** Convergence table for pdeg = ',num2str(pdeg),' ***']);
    disp('  1/numintervals |   ndofs   |    L2-error     |     H1-error   ')
    disp('-------------------------------------------------------------------')
    for i = 1:length(i_s)
        params.numintervals = 2*2^i_s(i);
        params.pdeg = pdeg;
        model = poisson_model(params);
        model_data = gen_model_data(model);
        sim_data = detailed_simulation(model, model_data);
        
        % project exact solution onto higher degree polynomial fem func
        p4_df_info = Fem.Lagrange.DefaultInfo(4, model_data.grid);
        uexact_h = interpol_global(p4_df_info, model.solution);
        uh = interpol_local(p4_df_info, @(grid, el, lcoord, params) sim_data.uh(el, lcoord, params));
        
        err = uh - uexact_h;
        l2errs(i) = l2_norm(err);
        h1errs(i) = h1_norm(err);
        ndofs(i) = sim_data.uh.ndofs;
        disp(['   ',num2str(1/params.numintervals,'%10.5e'),'   |    ',...
            num2str(ndofs(i),'%10.4d'),'   |   ',...
            num2str(l2errs(i),'%10.5e'),'   |   ',...
            num2str(h1errs(i),'%10.5e')]);
    end
end

disp('press enter to continue');
pause;

% Stokes: Mini-Element
fprintf('\n');
disp('(b) Part 1: Mini-FE for Stokes problem');
disp('======================================');

params = [];
params.mesh_number = 1;
model = laminar_flow_model(params);
model.has_nonlinearity = 0; % linear
model.df_type = 'stokes_mini'; % Mini Element

% generate and plot grid:
model_data = gen_model_data(model);
disp('plotting grid ...');
figure, plot(model_data.grid);
axis tight
axis equal
title('Grid for Stokes problem');

% simulation and plot:
sim_data = detailed_simulation(model, model_data);
disp('plotting solution ...');
plot_params = [];
plot_params.modes = {'velocity_abs'}; % for now, plot only velocity
plot_sim_data(model, model_data, sim_data, plot_params);

disp('press enter to continue');
pause;

% Taylor-Hood
fprintf('\n');
disp('(b) Part 2: Taylor-Hood FE for Stokes problem');
disp('=============================================');

% change model, new computations
model.df_type = 'taylor_hood'; % P2-P1
model_data = gen_model_data(model);
sim_data = detailed_simulation(model, model_data);
disp('plotting solution ...');
plot_sim_data(model, model_data, sim_data, plot_params);

disp('press enter to continue');
pause;
close all

% demonstrate precomputation of operator components (Fem.OperatorsDefault
% class)
fprintf('\n');
disp('Operator components can be precomputed (parameter decomposition):');

% finer mesh:
params.mesh_number = 4;
model = laminar_flow_model(params);
model.has_nonlinearity = 0;
model_data = gen_model_data(model);
sim_data = detailed_simulation(model, model_data);
disp(['Simulation time (full assembly): t = ', num2str(sim_data.total_time)]);

fprintf('Precompute operator components in ...');
tic;
model_data.operators.assemble_components(model, model_data);
fprintf([' t = ', num2str(toc), '\n']);

sim_data = detailed_simulation(model, model_data);
disp(['Simulation time (using stored data): t = ', num2str(sim_data.total_time)]);

% erase components:
model_data.operators.clear_components();

disp('press enter to continue');
pause;

% nonlinear model:
fprintf('\n');
disp('(c): Nonlinear Stokes problem (CFD benchmark)');
disp('=============================================');
params.mesh_number = 1;
model = laminar_flow_model(params);
model_data = gen_model_data(model);
model.verbose = 4;
disp(['Simulation for RE = ', num2str(model.get_reynolds_number(model))]);
sim_data = detailed_simulation(model, model_data);
disp('plotting solution ...');
plot_sim_data(model, model_data, sim_data, plot_params);

% higher reynolds number:
fprintf('\n');
model = set_mu(model, 2);
disp(['Simulation for RE = ', num2str(model.get_reynolds_number(model))]);
sim_data = detailed_simulation(model, model_data);
disp('plotting solution ...');
plot_sim_data(model, model_data, sim_data, plot_params);

disp('press enter to continue');
pause;

% finer mesh for demonstrating caching
fprintf('\n');
disp('Speed-up assembly by caching:');
disp('Perform simulation on finer mesh:');
fprintf('(without caching) ');
params.mesh_number = 4;
model = laminar_flow_model(params);
model = set_mu(model, 2);
model.verbose = 3; % less information
cache('limit', 2e8);
model_data = gen_model_data(model);
detailed_simulation(model, model_data);

fprintf('(with caching) ');
model.disable_caching = 0;
sim_data = detailed_simulation(model, model_data);

disp('press enter to continue');
pause;
close all

% plotting functionality
fprintf('\n');
disp('Show all plot modes for Stokes problems:');
plot_params.modes = {...
    'pressure', 'velocity', 'velocity_xcomp', 'velocity_ycomp', ...
    'velocity_abs', 'velocity_vec', 'velocity_vec_plus_pressure'};
plot_params.vector_distance = 0.03;
plot_sim_data(model, model_data, sim_data, plot_params);

disp('press enter to continue');
pause;

% compute benchmark values and compare to reference
fprintf('\n');
disp('Compute drag and lift coefficients on different grids:');
disp('Reference values (benchmark): C_D = 5.57953523384, C_L = 0.010618948146');
fprintf('\n');
disp('       h       |      ndofs     |     C_D    |     C_L');
disp('----------------------------------------------------------');

for i = 1:4
    params.mesh_number = i;
    model = laminar_flow_model(params);
    model_data = gen_model_data(model);
    model = set_mu(model, 1);
    model.fp_maxiter = 17;
    sim_data = detailed_simulation(model, model_data);
    [C_D, C_L] = stokes_compute_drag_and_lift(model, model_data, sim_data);
    disp(['   ' ,num2str(max(max(model_data.grid.EL)), '%.3e'), ...
        '   |    ', num2str(model_data.df_info.ndofs, '%.3e'),'   |   ', ...
        num2str(C_D, '%.4f'),'   |   ', num2str(C_L, '%.4f')]);
end

disp('press enter to terminate (close all windows)');
pause;
close all