buckley_leverett_model.m

function model = buckley_leverett_model(params)
% function model = buckley_leverett_model(params)
%
% optional fields of params:
%  model_size:   - 'big'    simulations are done on a larger grid with very
%                small time steps
%                - 'small'  simulations are done on a smaller grid with
%                bigger time steps
%  model_type:   - 'buckley_leverett' -- A Buckley-Leverett problem
%
  model = nonlin_evol_model_default;
  model.name = 'buckley_leverett_model';
  model.rb_problem_type = 'bl';
  %% data that is definitely used outside of the detailed simulations
  %  model.mu_names = {'diff_m', 'diff_p', 'hill_height'};
  %  model.mu_ranges = {[0,0.5],[0.01,100],[0,0.3]};
  model.mu_names = {'diff_k0'};
  model.mu_ranges = {[0.01, 0.5]};


  %% data that might be of interest outside of the detailed simulations
  model.T  = 0.2;
  model.nt = 80;

  %% finite volume settings
%  model.diffusivity_ptr            = @diffusivity_buckley_leverett;
%  model.diffusivity_derivative_ptr = @diffusivity_buckley_leverett_derivative;
  model.diffusivity_ptr            = @diffusivity_buckley_leverett_simple;
  model.diffusivity_derivative_ptr = @diffusivity_buckley_leverett_simple_derivative;

  model.num_conv_flux_ptr       = @fv_num_conv_flux_engquist_osher;
  model.num_diff_flux_ptr       = @fv_num_diff_flux_gradient;
  model.flux_linear             = 0;
  model.name_flux               = 'exponential';
  model.geometry_transformation = 'none';

  model.newton_solver        = 1;
  model.collect_newton_steps = false;
  model.newton_steps         = 60;


  model.operators_diff_implicit    = @fv_operators_diff_implicit_gradient;
  model.operators_conv_implicit    = @fv_operators_zero;
  model.operators_neumann_implicit = @fv_operators_zero;
  model.fv_expl_conv_weight  = 1.0;
  model.fv_impl_diff_weight  = 1.0;
  model.implicit_gradient_operators_algorithm = ...
    @fv_frechet_operators_diff_implicit_gradient;

  %    lxf_lambda               = 1.25;
  %    lxf_lambda               = fv_search_max_lxf_lambda([],    


  %% define the grid
  model = unitcube(model);
  model.xnumintervals = 100;
  model.ynumintervals = 100;

  % pointer to function in rbmatlab/datafunc/dirichlet_values
%  model.dirichlet_values_ptr = @dirichlet_values_uplow;
%  model.c_dir                = 0.2;
%  model.dir_box_xrange       = [-1 2];
%  model.dir_box_yrange       = 0.5 + [1/320 3/320];
%  model.c_dir_left           = 1;
%  model.c_dir_right          = 1;
%  model.c_dir_up             = 0;
%  model.c_dir_low            = 0.5;
%  model.dir_middle           = 0.7;
  model.dirichlet_values_ptr = @(glob, params) params.c_low * ones(length(glob),1);

  % model.dirichlet_values_ptr = @dirichlet_values_leftright;
  model.c_dir_left  = 0.0;
  model.c_dir_right = 0.0;
  model.dir_middle  = 0.5;
%  model.c_dir = 0.1;

  % pointer to function in rbmatlab/datafunc/neumann_values
  %model.neumann_values_ptr = @neumann_values_homogeneous;
  model.neumann_values_ptr = @neumann_values_zero;
  model.c_neu              = 0.00;

  % name of function in rbmatlab/datafunc/init_values/
%  model.init_values_ptr = @init_values_transformed_blobs;
%  model.init_values_ptr = @init_values_bars;
  model.init_values_ptr = @init_values_bl;
  % parameters for data functions

  model.filecache_velocity_matrixfile_extract = 0;

  % settings for CRB generation
  model.CRB_generation_mode = 'param-time-space-grid';
  model.RB_generation_mode  = 'greedy_refined';
  model.RB_refinement_mode  = 'uniform';
  model.RB_stop_max_refinement_level = 3;
  model.RB_val_rand_seed    = 12345;
  model.refinement_theta    = 0.6;
  model.Mmax    = 250;
  model.sMmax   = {50, 50};
  model.M       = 40;
  model.Mstrich = 5;
  model.ei_stop_on_Mmax = 1;
  % build ei-basis with WHOLE trajectory
  model.ei_target_error = 'linfty';
  model.RB_stop_max_val_train_ratio = 1.4;

  % target accuracy epsilon:

  % decide, whether estimator or true error is error-indicator for greedy
  model.RB_error_indicator = 'error'; % true error
  % RB_error_indicator = 'estimator'; % Delta from rb_simulation
  % RB_error_indicator = 'ei_estimator_test'; % Delta from rb_simulation testet
  % against true error

  model.divclean_mode    = false;
  model.flux_quad_degree = 1;

%  model.conv_flux_ptr             = @conv_flux_brooks_corey;
%  model.conv_flux_derivative_ptr  = @conv_flux_brooks_corey_derivative;
  model.conv_flux_ptr             = @conv_flux_brooks_corey_simple;
  model.conv_flux_derivative_ptr  = @conv_flux_brooks_corey_simple_derivative;

  % set all to dirichlet-boundary by specifying "rectangles", all
  % boundary within is set to boundary type by bnd_rect_index
  model.bnd_rect_corner1 = [-0.999999, 0; ...
                             0,       -0.9999999 ];
  model.bnd_rect_corner2 = [ 2,        0.999999; ...
                            0.999999, 2 ];
  % -1 means dirichlet, -2 means neumann
  model.bnd_rect_index   = [ -2, -2 ];

  if isfield(params, 'separate_CRBs') && params.separate_CRBs
    model.separate_CRBs = params.separate_CRBs;
    savepath_infix      = '';
  else
    savepath_infix      = '_one_CRB';
  end

  model.ei_detailed_savepath = [model.name,'_',params.model_type,...
                                savepath_infix,'_ei_data_interpol'];
  model.ei_operator_savepath = [model.name,'_',params.model_type,...
                                savepath_infix,'_ei_operators_interpol'];
  model.RB_detailed_val_savepath = [model.name,'_',params.model_type,...
                                savepath_infix,'_RB_val_detailed'];
  model.RB_detailed_train_savepath = [model.name,'_',params.model_type,...
                                savepath_infix,'_RB_train_detailed'];
%  model.RB_detailed_train_savepath = model.ei_detailed_savepath;

  if isfield(params, 'model_size')
    if strcmp(params.model_size, 'big') == 1
      model.xnumintervals = 200;
      model.ynumintervals = 200;
    end
    if strcmp(params.model_size, 'small') == 1
      model.xnumintervals = 50;
      model.ynumintervals = 50;
    end
  end

  if isfield(params, 'model_type')
    if strcmp(params.model_type, 'buckley_leverett')
      model.diff_K                = 1.0000;
      model.bl_k                  = 1/0.20;
      model.bl_mu1                = 5;
      model.bl_mu2                = 5;
      model.bl_lambda             = 0.3;
      model.newton_regularisation = false;
      model.c_low                 = 0.02;
      model.c_high                = 1.0;
      model.T                     = 0.3;
      model.mu_names              = {'diff_K','c_low','bl_lambda'};
      model.mu_ranges             = {[0.1,2],[0.0,0.04],[0.1,0.4]};
      model.adaptive_time_split_mode = false;
      model.ei_stop_epsilon_first = 5e-3;
      model.ei_stop_epsilon       = 1e-6;
      model.minimum_time_slice    = 5;
      model.time_split_Mmax       = 75;
      model.Mmax                  = Inf;
      model.extend_crb            = true;
      model.skip_search           = true;

%      model.mu_names              = {'bl_mu2','diff_K','bl_lambda'};
%      model.mu_ranges             = {[1,50],[0,2],[0.1,0.4]};
%      model.ei_numintervals       = [8,4,3,3];
%      model.RB_numintervals       = [8,4,3,3];
      model.ei_numintervals       = [4,3,3];
      model.RB_numintervals       = [2,1,1];
      model.nt                    = 60;
      model.RB_stop_Nmax          = 250;
      model.RB_stop_epsilon       = 5e-4;
      model.RB_stop_timeout       = 5*60*60;
      model.bnd_rect_index        = [-1,-1];
      model.separate_CRBs         = false;
      model.ei_space_operators    = { model.L_E_local_ptr , model.L_I_local_ptr };
      model.c_width = 0.2;
      model.c_height = 0.4;

%      if isfield(params, 'use_laplacian') && params.use_laplacian
%        model.laplacian_derivative_ptr = @(glob, U, model) ...
%          ones(size(U));
%        model.laplacian_ptr            = @(glob, U, model) ...
%          U;
%        model.diffusivity_ptr          = @(glob, U, model) ...
%          struct('K',1,'epsilon',0);
%        model.diffusivity_derivative   = @(glob, U, model) ...
%          struct('K',1,'epsilon',0);
%        model.implicit_gradient_operators_algorithm = ...
%          @fv_operators_diff_implicit_gradient;
%      end
%    elseif strcmp(params.model_type, 'eoc')
%      model.xnumintervals        = 10;
%      model.ynumintervals        = 10;
%      model.diff_k0              = 0.05;
%      model.diff_m               = 0;
%      model.nt                   = 50;
%      model.init_values_ptr      = @exact_function_heat_equation;
%      model.dirichlet_values_ptr = @exact_function_heat_equation;
%      model.bnd_rect_index       = [ -1, -1 ];
%    elseif strcmp(params.model_type, 'eoc_nonlinear')
%      model.xnumintervals        = 10;
%      model.ynumintervals        = 10;
%      model.newton_regularisation = 0;
%      model.diff_k0              = 0.0;
%      model.diff_m               = 1;
%      model.diff_p               = 1.0;
%      model.nt                   = 150;
%      model.T                    = 1.0;
%      model.init_values_ptr      = @exact_function_plaplace;
%      model.dirichlet_values_ptr = @exact_function_plaplace;
%      model.data_const_in_time   = false;
%      model.bnd_rect_index       = [ -1, -2 ];
%    elseif strcmp(params.model_type, 'buckley_leverett')
    end
  end

  model.model_type = params.model_type;

  model.implicit_nonlinear = false;

  model.data_const_in_time = true;

  model = model_default(model);

  if ~isfield(model, 'ei_space_operators')
    model.ei_space_operators = { model.L_E_local_ptr, model.L_I_local_ptr };
  end

  if isfield(params, 'verbose')
    model.verbose = params.verbose;
  else
    model.verbose = 0;
  end
  model.debug   = 1;

end
% vim: set et sw=2:
%| \docupdate