newton_oo_model.m

function model = newton_oo_model(params)
% function model = newton_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
%                - 'linear_diffusion' -- A linear parabolic problem with
%                    constant diffusion
%                `` \partial_{t} u - \nabla \cdot ( \mu_1 \nabla u ) = 0 ``
%                - 'exponential_diffusion' -- A non-linear parabolic problem
%                    with exponential function as diffusion coefficient
%                `` \partial_{t} u 
%                   - \nabla \cdot ( k_0 + \mu_1 u^{\mu_2} \nabla u ) ``
%                - 'eoc' -- A linear diffusion problem with known exact
%                   solution for validation of numerical scheme
%                - 'eoc_nonlinear' -- A nonlinear diffusion problem with known
%                   exact solution for validation of numerical scheme
%

  if nargin == 0
    params = [];
  end

  if ~isfield(params, 'model_size')
    params.model_size = 'normal';
  end

  if ~isfield(params, 'model_type')
    params.model_type = 'exponential_diffusion';
  end

  model = NonlinEvol.descr_default;
  model.name = 'newton_oo_model';
  model.rb_problem_type = 'NonlinEvol';

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

  %% time domain
  model.T  = 1.0;
  model.nt = 80;

  %% model size, affecting grid & time domain
  if isfield(params, 'model_size')
    if isequal(params.model_size, 'big')
      model.xnumintervals = 200;
      model.ynumintervals = 200;
    elseif isequal(params.model_size, 'small')
      model.xnumintervals = 50;
      model.ynumintervals = 50;
      model.nt = 25;
    end
  end

  %% boundary settings
  % 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 ];

  %% diffusivity
  model.diffusivity_ptr            = @diffusivity_exponential;
  model.diffusivity_derivative_ptr = @diffusivity_exponential_derivative;
  model.diff_k0                    = 0.00;
  model.diff_m                     = 1.0;
  model.diff_p                     = 2.0;

  %% fluxes
  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.k       = 0.005;
  model.gravity = 0;

  %% operators
  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  = 0.0;
  model.fv_impl_diff_weight  = 1.0;
  model.implicit_gradient_operators_algorithm = ...
    @fv_frechet_operators_diff_implicit_gradient;

  model.implicit_nonlinear = true;

  model.data_const_in_time = false;

  %% dirichlet values
  model.dirichlet_values_ptr = @dirichlet_values_leftright;
  model.c_dir_left  = 0.1;
  model.c_dir_right = 0.6;
  model.dir_middle  = 0.5;

  %% neumann values
  model.neumann_values_ptr = @neumann_values_homogeneous;
  %model.neumann_values_ptr = @neumann_values_outflow;
  model.c_neu              = 0.00;

  %% init values
  model.init_values_ptr = @init_values_bars;
  model.c_init = 0.1;

  %% velocity/convection
  model.filecache_velocity_matrixfile_extract = 0;
  model.divclean_mode    = false;
  model.flux_quad_degree = 1;

  model.conv_flux_ptr             = @conv_flux_linear;

  model.velocity_ptr              = @velocity_transport;
  model.transport_x               = 0.1;
  model.transport_y               = 0.1;

  model.conv_flux_derivative_ptr  = @(glob, U, params) ...
                                     params.velocity_ptr(glob, params);


  %% solver
  model.newton_solver     = 1;
  model.newton_steps      = 60;

  model.newton_regularisation = 0;

  model.stencil_mode       = 'edge';
  model.local_stencil_size = 1;

  %% parameters
  model.mu_names = {'diff_k0'};
  model.mu_ranges = {[0.01, 0.5]};

  % default rb values
  model.parspacesize = 8;
  stop_Nmax    = 20;
  stop_epsilon = 1e-4;
  stop_timeout = 24*60*60;

  if isfield(params, 'model_type')
    if strcmp(params.model_type, 'linear_diffusion')
      model.mu_names        = {'diff_k0'};
      model.mu_ranges       = {[0.01, 0.5]};
    elseif strcmp(params.model_type, 'exponential_diffusion')
      model.diff_k0         = 0.0000;
      model.diff_m          = 0.01;
      model.bnd_rect_index  = [-2,-2];
      model.c_init          = 0.0001;
      model.newton_regularisation = 0;
      model.diff_p          = 2.0;

      model.T               = 1;
      model.nt              = 80;

      model.mu_names        = {'diff_p','diff_m','c_init'};
      model.mu_ranges       = {[1 5],[0 0.01],[0 0.2]};

      if isfield(params, 'use_laplacian') && params.use_laplacian
        model.laplacian_derivative_ptr = @(glob, U, model) ...
          real((model.diff_m*U.^(model.diff_p-1)));
        model.laplacian_ptr            = @(glob, U, model) ...
          real((1/model.diff_p * model.diff_m * U.^(model.diff_p)));
        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

      model.parspacesize          = [4, 1, 2];

    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.T                    = 1.0;
      model.nt                   = 150;

      model.init_values_ptr      = @exact_function_plaplace;
      model.dirichlet_values_ptr = @exact_function_plaplace;

      model.bnd_rect_index       = [ -1, -2 ];

    elseif strcmp(params.model_type, 'buckley_leverett')
      error('buckley leverett problem type is not implemented yet.');
    end
  end

  model.enable_error_estimator = 1;
  model.Mstrich                = 1;

  % error norm for residual in error estimation
  model.error_norm = 'l2';

  %% finalize the model
  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   = 0;


end

% vim: set et sw=2: