advection_diffusion_fv_model.m

function model = advection_diffusion_fv_model(params);
%function model = advection_diffusion_fv_model(params);
%
% advection-diffusion model 
% two overlapping divergence free velocity fields, uniform
% diffusivity. weight of x-velocity-field and diffusivity
% are the parameters. 
% Pure explicit discretization of convection with fv
% discretization, implicit discretization of diffusion. 
% In contrast to the advection_fv_output_model, 
% the velocity field is constant in time (but parameter
% dependent). The boundary conditions are decaying in time.
% see rb_tutorial steps 11 - end for usage
% Model to be used for the Luminy book chapter 
%
% possible fields of params:
%      coarse_factor: coarsening factor between 1 and 8 
%                     time-steps, gridsize are scaled down with this. 
%

% B. Haasdonk 10.2.2014

model = [];
if (nargin ==1) & (isfield(params,'coarse_factor'))
  coarse_factor = params.coarse_factor;
else
  coarse_factor = 1;
end;

model = lin_evol_model_default;
model.name = 'advection_diffusion_model';

% RB information:
%model.mu_names           = {'vx_weight','k'};
model.mu_names           = {'mu1','mu2'};
model.mu_ranges          = {[0 1],[0,1]};

% specification of the time information
model.T             = 1;
model.nt            = 1024/coarse_factor;
%model.save_time_indices = 0:16/coarse_factor:model.nt;
model.save_time_indices = 0:4/coarse_factor:model.nt;
%disp('nt to be adjusted later!');
model.dt = model.T/model.nt;
model.verbose = 9;
model.debug = 0;
%model.debug = 1; %
%if model.debug
%  disp('model with debugging turned on.');
%end;
model.axis_equal = 1;
model.error_estimation = 1;

% grid information
%model.grid_initfile = 'rectangle_triagrid.mat';
model.gridtype = 'rectgrid';
model.xnumintervals = 256/coarse_factor;
model.ynumintervals = 128/coarse_factor;
model.xrange = [0,2];
model.yrange = [0,1];
% set completely dirichlet boundary on 4 edges
%model.bnd_rect_corner1 = [0,0;0,0;0,1;2,0]'-2*eps;
%model.bnd_rect_corner2 = [0,1;2,0;2,1;2,1]'+2*eps;
%model.bnd_rect_index   = [-1,-1,-1,-1];
model.bnd_rect_corner1 = [0,0]'-2*eps;
model.bnd_rect_corner2 = [2,1]'+2*eps;
model.bnd_rect_index   = [-1];

% set Dirichlet on upper and left, set outflow on bottom and right
%model.bnd_rect_corner1 = [0,0;0,0;0,1;2,0]'-2*eps;
%model.bnd_rect_corner2 = [0,1;2,0;2,1;2,1]'+2*eps;
%model.bnd_rect_index   = [-1,-2,-1,-2];

model.element_quadrature = @triaquadrature;
model.intersection_quadrature = @intervalquadrature;
model.dim_U = model.xnumintervals * model.ynumintervals * 2;

% Dirichlet and Initial data
model.dirichlet_values_ptr = @dirichlet_values_affine_decomposed;
model.dirichlet_values_coefficients_ptr = @my_dirichlet_values_coefficients;
model.dirichlet_values_components_ptr = @my_dirichlet_values_components;
% dummy one dirichlet_values:
%model.dirichlet_values_coefficients_ptr = @(params) 1;
%model.dirichlet_values_components_ptr = @(glob,params) {ones(1, ...
%                                                 size(glob,2))};
%model.neumann_values_ptr = @neumann_values_outflow;
model.init_values_ptr = @init_values_affine_decomposed;
model.init_values_coefficients_ptr = @my_dirichlet_values_coefficients_t0;
model.init_values_components_ptr = @my_dirichlet_values_components;
model.cone_number = 3; % number of components = cones 
model.cone_range = [0,1]; % x-range of cone-support 
model.cone_weight = 1; % leftmost cone with full weight
model.cone_amplitude = 1; % full amplitude
model.velocity_ptr = @velocity_affine_decomposed;
model.velocity_coefficients_ptr = @my_velocity_coefficients;
model.velocity_components_ptr = @my_velocity_components;
%model.vx_weight = 1.0;
%model.vx_weight = 0.75;
%model.vy_weight = 1.0;
model.vy_weight = 0.75;
model.conv_flux_ptr = @conv_flux_linear;
model.diffusivity_ptr = @diffusivity_homogeneous;
%model.k  = 1e-2;

model.rb_init_data_basis = @RB_init_data_basis;

% POD-Greedy settings:
model.RB_generation_mode = 'greedy_uniform_fixed';
model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
model.RB_stop_epsilon = 1e-9; 
model.RB_stop_timeout = 60*60; % 1 hour
model.RB_stop_Nmax = 50;
%model.RB_stop_max_val_train_ratio = inf;
model.filecache_ignore_fields_in_model = {'N','Nmax'};
model.filecache_ignore_fields_in_detailed_data = {'RB_info'};

% decide, whether estimator or true error is error-indicator for greedy
% params.RB_error_indicator = 'error'; % true error
model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
% choose a 5x5 uniform grid
model.RB_numintervals = 4 * ones(size(model.mu_names));

% test parameters
%model.RB_detailed_test_savepath = 'test_data_100';
%model.RB_test_size = 1000;

% perhaps these are redundant later...
model.divclean_mode = 0;
model.flux_quad_degree = 1;
model.flux_linear = 1;

model.init_values_algorithm = @disc_init_values;
model.init_values_qdeg = 0;
model.pdeg = 0; % FV schemes with piecewise constant ansatz functions
model.evaluate_basis = @fv_evaluate_basis;
model.l2project = @l2project;
model.local_mass_matrix = @fv_local_mass_matrix_tria; % triangular grid
model.mass_matrix = @fv_mass_matrix; 
model.ndofs_per_element = 1;
model.plot = @fv_plot; % then also plot_sequence will work
model.set_mu = @my_set_mu;

% for use of old datafunctions, use the following:

%model.operators_ptr = @fv_operators_implicit_explicit;
% for use of new data function access, use now
model.operators_conv_explicit = @fv_operators_conv_explicit_engquist_osher;
model.operators_conv_implicit = @fv_operators_zero;
model.operators_diff_explicit = @fv_operators_zero;
model.operators_neumann_implicit = @fv_operators_zero;
model.operators_neumann_explicit = @fv_operators_zero;
model.operators_diff_implicit = @fv_operators_diff_implicit_gradient;
%model.operators_neumann_explicit = @fv_operators_neumann_explicit;
%model.operators_output = @fv_operators_output;

%model.flux_linear = 1;
model.data_const_in_time = 0; % time varying data...
%model.diffusivity_ptr = ...;
%model.neumann_value_ptr = @...
model.affinely_decomposed = 1; % data functions allow affine
                               % parameter decomposition

model.compute_output_functional = 0; % turn off computation of output
%model.output_function_ptr = @output_function_box_mean;
%model.sbox_xmin = 1;
%model.sbox_xmin = 1;
%model.sbox_xmax = 2;
%model.sbox_ymin = 0;
%model.sbox_ymax = 0.5;
%model.sbox_ymax = 0.5;

%model.disc_element_mean_ptr = @fv_element_mean;

model = model.set_mu(model,[0,0]);

return; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% auxiliary functions used as pointers above:
%%%%% check later, if they are of general interest to be exported
%%%%% as standalone functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% new function syntax: global coordinates as rows in glob_coord,
% time and parameter variables in params

% dirichlet-data: convex combination of equidistant cones, decaying
% in time
function res = my_dirichlet_values_coefficients_t0(params);
params.t = 0;
res = my_dirichlet_values_coefficients(params);

function res = my_dirichlet_values_coefficients(params);
% res is a vector of coefficients
res = (1-params.t);

function res = my_dirichlet_values_components(glob,params);
% res is a cell-array of row-vectors of global evaluations
Q_0 = 3;
delta_cone = (params.cone_range(2)-params.cone_range(1))/(Q_0+1);
cone_pos_x = delta_cone * (1:Q_0)+ params.cone_range(1);
q = 2;
res = cell(1,1);
res{1} = 1-min(sqrt((glob(:,1)-cone_pos_x(q)).^2+...
            (glob(:,2)-1).^2)*(Q_0+1),1);
res{1} = res{1} * params.cone_amplitude;

% velocity field: overlap of y- and x parabolic profiles
function res = my_velocity_coefficients(params);
%res = [params.vx_weight, params.vy_weight]*(1-params.t);
res = [params.vx_weight, params.vy_weight]*0.5; %*(1-params.t);

function res = my_velocity_components(glob,params);
resx = [5*(1-glob(:,2).^2),...
        zeros(size(glob,1),1)];
%resy = [zeros(1,size(glob,2));...
%        -2*((1-glob(1,:).^2 .*(glob(1,:)>=0) & (glob(1,:)<=1)...
%            ))];
resy = [zeros(size(glob,1),1),...
        -1*((4-glob(:,1).^2))];
res = {resx, resy};

% realize coupling of mu1 => vweight, mu2*0.03 => k
function model = my_set_mu(varargin)
model = set_mu_default(varargin{:});
model.vx_weight = model.mu1;
model.k = model.mu2*0.03;