buckley_leverett.m

Go to the documentation of this file.

% small script demonstrating a buckley leverett problem
% discretization and RB model reduction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Select here, what is to be performed with the richards data %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step = 1 % single detailed simulation with given data and plot. Run
          % this with varying parameters mu until sure that scheme
          % is stable. Modify dt or the data-functions accordingly,
          % until a nice parameter-domain with uniformly stable
          % detailed scheme is obtained.
%step = 2 % generate colateral reduced basis of L_E operator
%step = 3 % compare single detailed simulation with and without
          % interpolated space operator
%step = 4 % generate dummy reduced basis from single trajectory and check, if
          % ei_interpolation with projection on this space maintains
          % result. A simple reduced simulation can also be
          % performed. All results should be visually identical
%step = 5 % generate reduced basis
%step = 6 % time measurement of reduced simulation and
          % use reduced basis in rb_demo_gui
%step = 7 % plot trajectories for some parameters mu at few time-instants
%step = 8 % do a coupled error plot

steps = [0,2,5,7,9]
%steps = [0,1]
%steps = [0,7.5]

warning('off', 'RBMatlab:rb_simulation:runtime');

imsavepath = fullfile(rbmatlabresult, 'images');

for si=1:length(steps)

step = steps(si);

params.model_type    = 'buckley_leverett';
%params.model_type     = 'eoc_nonlinear';
%params.model_size    = 'big';
%params.model_size    = 'medium';
params.model_size    = 'small';
%  params.model_type    = 'nonlinear';
params.separate_CRBs = false;

% test parameters
params.RB_detailed_test_savepath = 'test_data_100';
params.RB_test_size              = 100;

params.verbose = 1;
params.debug   = 1;
params.step7_outputfile = fullfile(rbmatlabresult, [params.model_type, '_step7_output']);
params.step8_outputfile = fullfile(rbmatlabresult, [params.model_type, '_step8_output']);
extrapar.Mstrich = 20;
extrapar.N       = 25;
extrapar.M       = 50;

% email to which status messages about processing mu vectors is sent
% params.email   = 'mdrohmann@uni-muenster.de';

%% parameters for visualization
plot_params.show_colorbar = 1;
plot_params.colorbar_mode = 'EastOutside';
plot_params.no_lines      = true;
plot_params.plot_type     = 'mesh';
% plot_params.clim          = [0, 0.5];
% plot_params.clim_tight    = true;
plot_params.plot          = @fv_plot;
plot_params.bind_to_model = true;
plot_params.yscale_uicontrols = 0.6;

% output-filenames in rbmatlabresult
CRBfname      = ['bl_', params.model_type, '_CRB_interpol.mat'];
detailedfname = ['bl_', params.model_type, '_detailed_interpol.mat'];

switch step
 case 0 % initialize model data
   model            = buckley_leverett_model(params);
   model.force_delete = true;
   model.debug = false;
   params.mu_ranges = model.mu_ranges;
   params.mu_names  = model.mu_names;

   model_data = gen_model_data(model);
 case 1 % single detailed simulation and plot
   step1_detailed_simulation;
 case 1.5
   tmppar.model_type         = 'eoc';

   model = newton_model(tmppar);
   model.debug = 1;
   model.newton_solver = 0;
   model.T = 1;
   model.dt = model.T/model.nt;
   model.verbose = 21;
   err = zeros(model.nt+1,7);
   eoc = zeros(4, 1);
   for i = 1:7
     model_data   = gen_model_data(model);

     sim_data     = detailed_simulation(model, model_data);

     ffe = @exact_function_heat_equation;

     exact_data.U = zeros(model_data.grid.nelements, model.nt+1);

     for tt = 1:model.nt+1
       model.t = (tt-1) * model.dt;
       exact_data.U(:,tt) = ffe([model_data.grid.CX(:), model_data.grid.CY(:)], model);
     end
     model.t = 0;

     err(:,i)       = fv_l2_error(sim_data.U, exact_data.U, model_data.W);
     if(i>1)
       eoc(i-1)   = log(max(err(:,i-1))/max(err(:,i)))/log(2);
       disp(['Error: ', num2str(err(end,i)), 'EOC: ', num2str(eoc(i-1))]);
       pause;
     end

     model.xnumintervals = round(model.xnumintervals.*2);
     model.ynumintervals = round(model.ynumintervals.*2);
   end
 case 1.75
   tmppar.model_type         = 'eoc_nonlinear';
   tmppar.local_mass_matrix  = @fv_local_mass_matrix_rect;
   tmppar.element_quadrature = @cubequadrature;
   tmppar.evaluate_basis     = @fv_evaluate_basis;
   tmppar.pdeg = 0;

   model = newton_model(tmppar);
   %model.debug = 1;
   % model.newton_regularisation = 0;
   %model.verbose         = 21;
   model.dt              = model.T/model.nt;
   model.ynumintervals = 1;
   err = zeros(model.nt+1,5);
   eoc = zeros(4, 1);
   for i = 1:6
     disp(['Computing EOC step no. ', num2str(i)]);
     model_data   = gen_model_data(model);

     sim_data     = detailed_simulation(model, model_data);

     ffe = @exact_function_plaplace;

     width = model_data.grid.X(2) - model_data.grid.X(1);

     Fe = @(einds, loc, grid, params) ffe([grid.X(einds) + loc(:,1)*width,...
                                           grid.Y(einds) + loc(:,2)*width], ...
                                           params);
     tmppar.diff_k0 = model.diff_k0;
     exact_data.U = zeros(model_data.grid.nelements, model.nt+1);
     for tt = 1:model.nt+1
       model.t = (tt-1) * model.dt;
       exact_data.U(:,tt) = ffe([model_data.grid.CX(:), ...
                                 model_data.grid.CY(:)], model);
%       exact_data.U(:,tt) = l2project(Fe, 2, model_data.grid, tmppar);
     end
     model.t = 0;

     err(:,i) = fv_l2_error(sim_data.U, exact_data.U, model_data.W);
     if(i>1)
       eoc(i-1)   = log(max(err(:,i-1))/max(err(:,i)))/log(2);
       disp(['Error: ', num2str(err(end,i)), ' EOC: ', num2str(eoc(i-1))]);
     end

     model.xnumintervals = model.xnumintervals.*2;
     model.nt = model.nt;
     model.dt = model.T / model.nt;
     %model.ynumintervals = model.ynumintervals.*2;
   end

 case 2 % empirical interpolation of space operator
   step2_empirical_interpolation;
 case 3 % compare detailed simulation with and without interpolation
   model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
   step3_detailed_ei_simulation;
 case 4 % construct dummy reduced basis by single trajectory and simulate
   model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
%   model = model.set_mu(model, [0.1,0.01, 0.4]);
   model = model.set_mu(model, [0.11,0.00, 0.12]);
   step4_dummy_reduced_basis;
 case 5 % reduced basis
%   load(fullfile(rbmatlabresult,CRBfname));
   step5_rb_generation;
 case 6
%   load(fullfile(rbmatlabresult,detailedfname));
   step6_demo_rb_gui;

 case 7
  tmp = load(fullfile(rbmatlabresult, detailedfname));
  detailed_data = tmp.detailed_data;
  mu_set = {[2,0.1,0.4]};
  clim   = [0.0, 1.0];
  plot_params.plot_type = 'mesh';
  plot_params.line_width = 1;
  timeinstants = [0,0.3,1.0]*model.T;
  boxed_snapshots = true;
  step10_plot_trajectories;
case 8
  tmp=load(fullfile(rbmatlabresult, detailedfname));
  detailed_data = tmp.detailed_data;

  model.Mstrich   = 0;
  % maximum number of reduced basis functions
  model.Nmax      = size(detailed_data.RB,2);
  % maximum number of collateral reduced basis functions
  model.Mmax      = cellfun(@(x)size(x,2), detailed_data.BM);
  % number of reduced basis sizes to be tested.
  Nsize           = 7;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 10;
  % sample of values for the coupling variable 'c' for which errors are
  % computed.
  csample     = [1];
  % axis description for plots with coupling variable 'c'
  cdescr      = '(N,M) = c * (Nmax, Mmax)';
  % size of mu vector test set.
  mu_set_size = 20;

  error_plots = { 'error_coupled' };

  model.M_by_N_ratio = 3;

  step7_error_landscape;
case 9
  columntitles = {'minslicesize', 'noofslices', 'averageM', 'minM' ,'maxM', 'averagedtime'};
  values       = zeros(6,1);
  values(1) = model.minimum_time_slice;
  values(2) = max(detailed_data.time_split_map(:,2));
  Msizes = cellfun(@(X) size(X,2), detailed_data.BM);
  values(3) = sum(Msizes) / values(2);
  values(4) = min(Msizes);
  values(5) = max(Msizes);
  values(6) = mean(output{1}.rtimes);
  print_datatable(fullfile(rbmatlabresult,'tableline'),columntitles,values)
otherwise
  error('step-number is unknown!');
end

end

%| \docupdate