1.13.10/doc/conservation__test_8m_source.html

% test of conservativity of galerkin projection by the linear-fv approach.

% motivation of a conservative RB scheme is demonstrated

% with this example: The reduced scheme is not conservative, the

% missing conservativity correlates with the l1-error. So

% a conservative RB-scheme may result in a more accurate RB-scheme


% Bernard Haasdonk 25.1.2008


% parameters defining a linear problem with zero boundary conditions

% (then mass conservation follows)


params = [];

% rectangular or triangular

params.gridtype = 'triagrid';

params.grid_initfile = '2dsquaretriang';


%params.gridtype = 'rectgrid';

%params.xnumintervals = 20;

%params.ynumintervals = 20;

%params.xrange = [0,1];

%params.yrange = [0,1];


% set all to neuman-boundary by specifying "rectangles", all

% boundary within is set to boundary type

params.bnd_rect_corner1 = [-1; ...

                            -1];

params.bnd_rect_corner2 = [ 2; ...

                            2];


% -1 means dirichlet, -2 means neumann

params.bnd_rect_index =   -2; % first is dirichlet, second neuman


% time discretization

%params.T = 1.0;

%params.nt = 200; % guess and adjust later


params.T = 1.0;

params.nt = 80; % guess and adjust later


% output settings

params.verbose = 10;

params.axis_equal = 1;

params.clim = [-1,1];


% data functions

% name of function in RBmatlab/datafunc/init_values

params.name_init_values = 'gradient_box';

% parameters for data functions

params.c_init_up = 1.0;

%params.c_init_lo = 1.0;

params.c_init_lo = 0.0;

%params.c_init_lo = -1.0;


% name of function in RBmatlab/datafunc/dirichlet_values

params.name_dirichlet_values = 'zero';

%params.c_dir = 0;

%params.name_dirichlet_values = 'homogeneous';

%params.c_dir = 0;


% name of function in RBmatlab/datafunc/neuman_values

%params.name_neuman_values = 'homogeneous';

%params.c_neu = 0;

params.name_neuman_values = 'zero';


%params.name_neuman_values = 'homogeneous';

%params.c_neu = 0;


% convective flux


params.name_flux = 'parabola';

% maximum velocity:

params.c = 1;

% for simplifying computation in case of linear flux:

params.flux_linear = 1;

params.flux_quad_degree = 1;


%params.name_flux = 'burgers_parabola';

% for simplifying computation in case of linear flux:

%params.flux_linear = 0;

%params.flux_quad_degree = 2;


% parameter of chosen conv-flux

%params.vmax = 0.30;

%params.vrot_angle = 0;

%params.vrot_angle = -pi/4;

%params.vrot_angle = -pi/10;

%params.flux_pdeg = 2; % burgers exponent


% rb-formulation

params.mu_names = {'c_init_lo','c_init_up'};

params.mu_ranges = {[-1,1],[-1, 1]};

%params.rb_problem_type = 'nonlin_evol';

params.rb_problem_type = 'lin_evol';


% finite volume settings

params.name_diffusive_num_flux = 'none';

params.name_convective_num_flux = 'lax-friedrichs';


% determine maximum lxf

params.lxf_lambda = 3.6

if ~isfield(params,'lxf_lambda')

  disp('Lxf-Lambda is not set!! Computing empirically maximum possible...');

  params.lxf_lambda = 0.9 * fv_search_max_lxf_lambda([],params);

  disp(['Found maximum lxF-lambda =', num2str(params.lxf_lambda)]);

  disp('Please set this value in your script for skipping search next time!');

  keyboard;

end;

% later set explicit:

%params.lxf_lambda = 1.25;

%params.name_convective_num_flux = 'enquist-osher';


% projection of analytical initial data to discrete function

params.init_values_algorithm = 'fv_init_values';

% get mass matrix for inner product computation from DOF-vectors

params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';

params.data_const_in_time = 1;

params.L_I_inv_norm_bound = 1;

params.L_E_norm_bound = 1;


% settings for operators (only lin_evol)

params.operators_algorithm = 'fv_operators_implicit_explicit';


% settings for empirical interpolation (only for nonlin_evol)

%params.implicit_operators_algorithm = 'fv_operators_implicit';

%params.L_E_local_name = 'fv_conv_explicit_space';

%par.range = params.mu_ranges;

%params.ei_numintervals = [4,4]; % 8x8 parameter grid

%params.Mmax = 150;

%params.ei_detailed_savepath = 'burgers_fv_5x5_detailed';

%params.ei_operator_savepath = 'burgers_fv_5x5_ei_operator';

%params.ei_time_indices = 1:params.nt;

%params.CRB_generation_mode = 'param-time-space-grid';

%params.ei_target_error = 'interpol';


% settings for reduced basis generation

params.detailed_simulation_algorithm = 'detailed_simulation';

%params.init_values_algorithm = 'fv_init_values';

%params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';

params.error_algorithm = 'fv_error';

%params.lxf_lambda = 1.0194e+003;

%params.data_const_in_time = 1;

params.error_norm = 'l2';

params.RB_extension_algorithm = 'RB_extension_PCA_fixspace';

%                   'RB_extension_max_error_snapshot'};

params.RB_stop_timeout = 2*60*60; % 2 hours

params.RB_stop_epsilon = 1e-5;

params.RB_error_indicator = 'error'; % true error

%params.RB_error_indicator = 'estimator'; % Delta from rb_simulation

params.RB_stop_Nmax = 100;

%params.RB_generation_mode = 'uniform_fixed';

params.RB_numintervals = 5;

params.RB_detailed_train_savepath = 'conservation_test';


%disp('detailed simulation:');


detailed_data.grid = construct_grid(params);

U = detailed_simulation([],params);

params.title = 'detailed simulation';

plot_element_data_sequence(detailed_data.grid, U, params);

M = inner_product_matrix(detailed_data.grid,params);

% 1-norm

norm1U = sum( abs(M * U) );


% reduced simulation based on single trajectory


params.RB_generation_mode = 'PCA_trajectory';

detailed_data = rb_detailed_prep(params);


params.N = ceil(size(detailed_data.RB,2)/8);


offline_data = rb_offline_prep(detailed_data,params);

reduced_data = rb_online_prep(offline_data,params);

simulation_data = rb_simulation(reduced_data,params);

Urb = rb_reconstruction(detailed_data, simulation_data);


params.title = 'reduced simulation';

plot_element_data_sequence(detailed_data.grid, U, params);

norm1Urb = sum( abs(M * Urb) );


figure;

plot([norm1U; norm1Urb]')

legend({'detailed','reduced'});

title('1-norm of detailed and approximate solution');

xlabel('time-index k');

ylabel('1-norm');


figure;

err = fv_l2_error(U,Urb,detailed_data.W);

plot(err);

title('l2-error between detailed and approximate solution');

xlabel('time-index k');

ylabel('l2-error');


% TO BE ADJUSTED TO NEW SYNTAX

%| \docupdate