newton.m

Go to the documentation of this file.

function [descr, rdescr, dmodel, rmodel] = newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
% function [descr, rdescr, dmodel, rmodel] = newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
% 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

%steps = [0,2,3,4,5,6]
if nargin == 0
  steps = [0,5,7,8];
end

if nargin <= 1
  combined           = true;
  M_by_N_ratio       = 0;
  noof_ei_extensions = 1;
  use_laplacian      = 1;
  model_size         = 'normal';
  random             = false;
  num_cpus           = 4;
end

if nargin <= 8
  Mstrich = 1;
end

warning('off', 'RBMatlab:rb_simulation:runtime');
warning('off', 'MATLAB:cholinc:ArgInfToBeRemoved');

imsavepath = fullfile(rbmatlabresult, 'images');

rmodel = [];

for si=1:length(steps)

step = steps(si);

params.use_laplacian = use_laplacian;
%params.model_type    = 'buckley_leverett';
params.model_type    = 'exponential_diffusion';
extrapar.Mstrich     = 20;

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

params.step7_outputfile = fullfile(rbmatlabresult, ['newton_', params.model_type, '_step7.mat']);
params.step8_outputfile = fullfile(rbmatlabresult, ['newton_', params.model_type, '_step8.mat']);

% test parameters

params.verbose = 6;
params.debug   = 0;

% 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     = 'contour';
% 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
detailedfname = ['newton_', params.model_type, '_', num2str(Mstrich), '_detailed_interpol.mat'];

switch step
 case 0 % initialize model data
   if ~isempty(getenv('MACHINEFILE'))
     machinefile = getenv('MACHINEFILE');
     FakeMPI.MPI_Init(machinefile, fullfile(rbmatlabhome, 'fakeMPI', 'run_parCompErrs.sh'));
     FakeMPI.MPI_Run;
   end

   descr                  = newton_oo_model(params);
   rdescr.rb_problem_type = descr.rb_problem_type;
   if random
     rdescr.train_sample_mode        = 'random';
     rdescr.train_num_intervals      = 800;
     rdescr.stop_max_val_train_ratio = inf;
   end

   rdescr.stop_timeout                          = Inf;
   rdescr.stop_Nmax                             = 300;
   rdescr.stop_Mmax                             = 600;
   rdescr.stop_epsilon                          = 1e-4;
   rdescr.val_size                              = 20;
   rdescr.stop_max_val_train_ratio              = 1.10;
   rdescr.ei_minimum_residual                   = 1e-10;
   rdescr.extension_M_by_N_r                    = M_by_N_ratio;
   rdescr.noof_ei_extensions                    = noof_ei_extensions;
   rdescr.maximum_temporary_error_growth_factor = 1e3;
   rdescr.train_num_intervals                   = [4,1,2];
   rdescr.max_refinement_level                  = 8;
   rdescr.Mmax_small                            = 30;
   rdescr.ei_epsilon_small                      = 1e-3;
   if combined
     if Mstrich >= 1
       rdescr.stop_Nmax  = 111;
       rdescr.Mmax_small = rdescr.Mmax_small + Mstrich;
     end
   else
     rdescr.train_num_intervals                 = [7,2,4];
     rdescr.Mmax_small = 0;
     if Mstrich == 0
       rdescr.indicator_mode = 'error';
       rdescr.stop_Nmax      = 100;
       rdescr.stop_Mmax      = 426;
       rdescr.stop_epsilon   = 0;
     end
     if Mstrich > 1
       rdescr.stop_Nmax      = 100;
     end
   end


   rdescr.refinement_theta   = 0.10;

   [dmodel, rmodel] = gen_models(descr, rdescr);

   rmodel.Mstrich   = Mstrich;
   if Mstrich == 0
     rmodel.enable_error_estimator = false;
   end

   dmodel.num_cpus   = num_cpus;
   rmodel.num_cpus  = num_cpus;

   params.mu_ranges = rmodel.mu_ranges;
   params.mu_names  = rmodel.mu_names;

   model_data = gen_model_data(dmodel);
 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.grid, model);
     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.grid, model);
     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;
   step4_dummy_reduced_basis;
 case 5 % reduced basis
%   load(fullfile(rbmatlabresult,CRBfname));

   detailed_data = gen_detailed_data(rmodel, model_data);
   save(fullfile(rbmatlabresult, detailedfname), 'detailed_data', 'rmodel');
 case 6
   tmp=load(fullfile(rbmatlabresult,detailedfname));
   detailed_data=tmp.detailed_data;

   step6_demo_rb_gui;
 case 7 % training-error landscape & time comparisons
  load(fullfile(rbmatlabresult,detailedfname));

  basis_gen.Mstrich   = 0;
  % 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     = 0.1:0.1:1;

  M_by_N_ratio = 1.0;
  % axis description for plots with coupling variable 'c'
  cdescr      = ['(N,M) = c * (Nmax, ',num2str(M_by_N_ratio), '*Mmax)'];
  % size of mu vector test set.
  mu_set_size = 20;

  error_plots = { 'train_set_coupled', ...
                  'error_coupled' };


  step7_error_landscape;

case 7.5
  load(params.step7_outputfile);
  for i = 1:length(output)
    plot_error_landscape(output{i});
  end
case 8
  disp('warning: takes a few hours!');
  load(fullfile(rbmatlabresult,detailedfname));

  M_by_N_ratio = 1.0;

  % number of reduced basis sizes to be tested.
  Nsize           = 7;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 12;
  % Mstrich values to be tested
  Mstrich_samples = [1,2,3,5,8,10,15,20];
  % sample of values for the coupling variable 'c' for which
  % estimators are computed.
  csample         = 0.1:0.1:1;
  % size of mu vector test set
  mu_set_size     = 30;
  % axis description for plots with coupling variable 'c'
  cdescr        = '(N,M) = c * (Nmax, Mmax)';

  step8_estimator_landscape;
case 8.5
  load(params.step8_outputfile);
  for i = 1:length(output)
    output{i}.bound = 1500;
    output{i}.stab_limit = 10000;
    plot_error_landscape(output{i});
  end
case 9
  load(fullfile(rbmatlabresult,detailedfname));
  model_data = detailed_data;
  testdir = 'richards_test_10';
  dirnames = testdir;
  model.Nmax = size(detailed_data.RB,2);

  % generate test data if not existent
  rand('state',123456);
  M = rand_uniform(10,params.mu_ranges);
  save_detailed_simulations(model,model_data,M,testdir)

  % load demo_nonlin_evol_detailed_data_LE_on_RB;
  reduced_data = gen_reduced_data(model,detailed_data);

  settings = load(fullfile(rbmatlabtemp,...
                  testdir,'settings.mat'));

  Msamples = 12;
  Nsamples = 7;
  % storage for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
  errs   = zeros(Msamples,Nsamples);
  model.RB_error_indicator='ei_estimator_test';

  Mstrich_set = 1:10:100;%[1,2,3,4,5,10,15,20];

  deltas = zeros(size(M,2),length(Mstrich_set));
  for mu_ind = 1:size(M,2);
    disp(['processing parameter vector ',num2str(mu_ind),...
          '/',num2str(size(M,2))]);
    for msi = 1:length(Mstrich_set)

      model = model.set_mu(model,settings.M(:,mu_ind));
      sim_data = load_detailed_simulation(mu_ind,settings.savepath,model);
      U_H = sim_data.U;

      mtemp = model;

      mtemp.M       = extrapar.M;
      mtemp.Mstrich = Mstrich_set(msi);
      mtemp.N       = extrapar.N;
      simulation_data = rb_simulation(mtemp,reduced_data);
      disp('.');

      deltas(mu_ind,msi) = simulation_data.Delta(end);
      simulation_data = rb_reconstruction(mtemp,detailed_data,simulation_data);
      error =  model.error_algorithm(simulation_data.U, sim_data.U, ...
                                     detailed_data.grid, mtemp);
      disp(['true error: ', num2str(max(error))]);

%      plot_sim_data(model,detailed_data,simulation_data,plot_params);

    end
    plot(Mstrich_set,deltas(mu_ind,:));
    keyboard;
  end

case 9.8
%  tmp = load(fullfile(rbmatlabresult, 'exponential_laplace', 'newton_exponential_diffusion_detailed_interpol_new'));
  tmp = load(fullfile(rbmatlabresult, detailedfname));
  detailed_data = tmp.detailed_data;
  mu_set = {[1,0.01,0.2],[1,0.01,0.0],[4,0.01,0.2],[4,0.01,0.0]};
  clim   = [0.0, 0.6];
  plot_params.plot_type = 'contour';
  plot_params.contour_lines = 7;
  plot_params.line_width = 2;
  timeinstants = [0.1,1.0]*model.T;
  boxed_snapshots = true;
  step10_plot_trajectories;
case 10
  %tmp = load(fullfile(rbmatlabresult,'strange.mat'));
  %Cmap = tmp.Cmap;
  tmp                  = load(fullfile(rbmatlabresult,detailedfname));
  detailed_data        = tmp.detailed_data;
  model                = tmp.model;
  plot_params.no_lines = 1;

  plot_params.no_axes                = 1;
  plot_params.transparent_background = 1;
  plot_params.show_colorbar          = 0;
  %plot_params.shrink_factor         = 1.13;
  %plot_params.colormap               = Cmap;
  plot_params.clim                   = [0.1,0.6];

  reduced_data     = gen_reduced_data(model, detailed_data);
  model.N = size(detailed_data.RB,2);
  model.M = size(detailed_data.QM,2);
  reduced_data = extract_reduced_data_subset(model, reduced_data);

  cparams.range = cell(length(model.mu_ranges)+1,1);
  cparams.range{1}  = [1, model.nt];
  cparams.range(2:end) = model.mu_ranges(:);
  cparams.numintervals = 2*ones(length(model.mu_ranges)+1,1);
  cgrid = cubegrid(cparams);

  parammat = cgrid.vertex;
  parammat(:,1) = round(parammat(:,1));


  for i=1:length(parammat)
    tstep = parammat(i,1);
    newmodel = model.set_mu(model, parammat(i,2:end));
    rb_sim_data = rb_simulation(newmodel, reduced_data);
    rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);

    % transformed version
    plot_element_data(model_data.grid, rb_sim_data.U(:,tstep), plot_params);
    Uline = rb_sim_data.U(subline,tstep);
    mustring = strrep(strrep(strrep(strrep(mat2str(parammat(i,2:end)),' ', '_'),'[',''),']',''),'.','p');

    fp = fullfile(imsavepath, params.model_type);
    if ~exist(fp,'dir')
      mkdir(fp);
    end
    fp = fullfile(fp, ['t_', num2str(tstep)]);
    if ~exist(fp,'dir')
      mkdir(fp);
    end
    fn = ['sample_mu_', mustring];
    disp(['Processing ', fn, ' ...']);

    tikzparams.filename      = fn;
    tikzparams.filepath      = fp;
    tikzparams.save_colorbar = 1;

    plot_as_tikzfile(model, tikzparams);

    close(gcf);
  end

  datfn = 'sublines';
  print_datatable(fullfile(fp, datafn), subline_title, subline_values);

case 11
  load('richards_affine_detailed_interpol');

  rdata = gen_reduced_data(model, detailed_data);

  model.N = size(detailed_data,2);
  model.M = 10;
  model.Mstrich = 1;
  rdata2 = extract_reduced_data_subset(model, rdata);
  model.Mstrich = 10;
  rdata1 = extract_reduced_data_subset(model, rdata);

  TM1 = rdata1.TM_local{1}(1:10);
  TM2 = rdata2.TM_local{1}(1:10);

  gl1 = rdata1.grid_local_ext{1};
  gl2 = rdata2.grid_local_ext{1};

  if(any(gl1.CX(TM1) ~= gl2.CX(TM2)))
    disp('CX differs');
  end
  if(any(gl1.CY(TM1) ~= gl2.CY(TM2)))
    disp('CY differs');
  end
  if(any(any(gl1.X(gl1.VI(TM1,:)) ~= gl2.X(gl2.VI(TM2,:)))))
    disp('X differs');
  end
  if(any(any(gl1.Y(gl1.VI(TM1,:)) ~= gl2.Y(gl2.VI(TM2,:)))))
    disp('Y differs');
  end
  if(any(any(gl1.SX(TM1,:) ~= gl2.SX(TM2,:))))
    disp('SX differs');
  end
  if(any(any(gl1.SY(TM1,:) ~= gl2.SY(TM2,:))))
    disp('SY differs');
  end
    if(any(gl1.EL(TM1,:) ~= gl2.EL(TM2,:)))
      disp('EL differs');
    end
    if(any(gl1.DC(TM1,:) ~= gl2.DC(TM2,:)))
      disp('DC differs');
    end
    if(any(gl1.NX(TM1,:) ~= gl2.NX(TM2,:)))
      disp('NX differs');
    end
    if(any(gl1.NY(TM1,:) ~= gl2.NY(TM2,:)))
      disp('NY differs');
    end
    if(any(gl1.ECX(TM1,:) ~= gl2.ECX(TM2,:)))
      disp('ECX differs');
    end
    if(any(gl1.ECY(TM1,:) ~= gl2.ECY(TM2,:)))
      disp('ECY differs');
    end
  %% neighbours
  NBI1 = gl1.NBI(TM1,:);
  NBI2 = gl2.NBI(TM2,:);
  if(any(gl1.CX(NBI1) ~= gl2.CX(NBI2)))
    disp('CX differs');
  end
  if(any(gl1.CY(NBI1) ~= gl2.CY(NBI2)))
    disp('CY differs');
  end
  if(any(any(gl1.X(gl1.VI(NBI1,:)) ~= gl2.X(gl2.VI(NBI2,:)))))
    disp('X differs');
  end
  if(any(any(gl1.Y(gl1.VI(NBI1,:)) ~= gl2.Y(gl2.VI(NBI2,:)))))
    disp('Y differs');
  end
  if(any(any(gl1.SX(NBI1,:) ~= gl2.SX(NBI2,:))))
    disp('SX differs');
  end
  if(any(any(gl1.SY(NBI1,:) ~= gl2.SY(NBI2,:))))
    disp('SY differs');
  end
  for i = 1:length(TM1)
  end

case 12 % plot some reduced basis functions
  tmp                  = load(fullfile(rbmatlabresult,detailedfname));
  detailed_data        = tmp.detailed_data;
  model                = tmp.model;

  plot_params.no_lines = 1;

  plot_params.axis_equal = true;
  plot_params.plot_type = 'patch';
  plot_params.no_axes                = 1;
  plot_params.transparent_background = 1;
  plot_params.show_colorbar          = 0;

  subline = find(detailed_data.grid.CY > 0.6 ...
                 & detailed_data.grid.CY < 0.6+detailed_data.grid.hmin);
  subline_title = {'X'};
  subline_values = grid.CX(subline)';

  for i=1:10
    plot_element_data(detailed_data.grid, detailed_data.RB(:,i), plot_params);
    Uline = detailed_data.RB(subline,i);
    subline_title = [ subline_title, ['rb',num2str(i) ] ];
    subline_values = [ subline_values; reshape(Uline, 1, length(Uline)) ];

    fp = fullfile(imsavepath, params.model_type);
    if ~exist(fp, 'dir')
      mkdir(fp);
    end
    fn = ['rb_function_', num2str(i)];
    disp(['Processing ', fn, ' ...']);

    tikzparams.filename      = fn;
    tikzparams.filepath      = fp;
    tikzparams.save_colorbar = 1;
    tikzparams.width         = 3.5;
    tikzparams.scaleaxis     = '';

    plot_as_tikzfile(model, tikzparams);

    close(gcf);
  end
  datafn = 'rb_sublines.dat';
  print_datatable(fullfile(fp, datafn), subline_title, subline_values);


otherwise
  error('step-number is unknown!');
end

end

%| \docupdate