richards_fv.m

Go to the documentation of this file.

% small script demonstrating the richards equation with explicit fv
% discretization and RB model reduction
% Bernard Haasdonk 14.8.2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 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 % generate error-landscape over varying N and M
          % can take several hours!!!
%step = 8 % do runtime comparisons between detailed and reduced simulations
%step = 10 % write out some nice pictures

steps = {0,1}
%steps = {0,1}

imsavepath = '/data/dune_work/private/m_droh01/images';

for si=1:length(steps)

step = steps{si};

%  params.model_type    = 'implicit_nonaffine_linear';
params.model_type    = 'linear_heat_trapezoidal';
%params.model_size    = 'big';
params.model_size    = 'small';
% 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 = 5;
params.debug   = 0;
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.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      = ['richards_fv_', params.model_type, '_CRB_interpol.mat'];
detailedfname = ['richards_fv_', params.model_type, '_detailed_interpol.mat'];

switch step
 case 0 % initialize model data
   model            = richards_fv_model(params);
   params.mu_ranges = model.mu_ranges;
   params.mu_names  = model.mu_names;
   mu_default = cellfun(@mean, model.mu_ranges);

   model_data = gen_model_data(model);
 case 1 % single detailed simulation and plot
   dmodel = model;
   step1_detailed_simulation;
 case 2 % empirical interpolation of space operator
   step2_empirical_interpolation;
 case 3 % compare detailed simulation with and without interpolation
   load(fullfile(rbmatlabresult,CRBfname));
   step3_detailed_ei_simulation;
 case 4 % construct dummy reduced basis by single trajectory and simulate
   tmp=load(fullfile(rbmatlabresult,CRBfname));
   detailed_data=tmp.detailed_data;
   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 % training-error landscape & time comparisons
  load(fullfile(rbmatlabresult,detailedfname));

  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      = size(detailed_data.BM{1},2);
  % number of reduced basis sizes to be tested.
  Nsize           = 7;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 12;
  % sample of values for the coupling variable 'c' for which errors are
  % computed.
  csample     = 0.1:0.1: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;

  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));

  % maximum number of reduced basis vectors
  model.Nmax = model.RB_stop_Nmax;
  % maximum number of collateral reduced basis vectors used for
  model.Mmax = size( detailed_data.BM{1}, 2 ) - extrapar.Mstrich;

  % 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:extrapar.Mstrich;
  % sample of values for the coupling variable 'c' for which
  % estimators are computed.
  csamples        = 0.1:0.1:1;
  % size of mu vector test set
  mu_set_size     = 20;
  % 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.W, 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 10
  load('datafiles/strange.mat');

  tmp                  = load(fullfile(rbmatlabtemp,detailedfname));
  detailed_data        = tmp.detailed_data;
  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.0,0.35];

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

  model.hill_height = 0;

  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 = newmodel.set_mu(model, parammat(i,2:end));
    rb_sim_data = rb_simulation(newmodel, reduced_data, params);
    rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);

    % transformed version
    plot_element_data(newmodel, model_data.grid, rb_sim_data.U(:,tstep), plot_params);

    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, ' ...']);

    set(gcf,'Color','none');
    set(gca,'Color','none');
    set(gcf,'InvertHardCopy','off');
    showaxes('hide');
    ranges=cell(2,1);
    ranges{1} = get(gca,'XLim');
    ranges{2} = get(gca,'YLim');

    export_fig(fullfile(fp, fn), '-pdf','-eps','-png');
    im = export_fig;

    colorbarfn = fullfile(imsavepath, params.model_type, 'colorbar');
    if ~exist([colorbarfn, '.pdf'], 'file')
      colorh = colorbar;
      yticks = get(colorh, 'YTick');
      set(colorh, 'YTick',[]);
      set(colorh, 'YTickLabel',[]);

      export_fig(colorbarfn, '-pdf', '-eps', '-png', colorh);
      cbim = export_fig(colorh);
      size(cbim)
    end

    tikzfile(newmodel, fp, fn, ranges, size(im), 5);
    close(gcf);

    % untransformed version
    fpu = fullfile(imsavepath, params.model_type, 'untransformed');
    if ~exist(fpu,'dir')
      mkdir(fpu);
    end
    fpu = fullfile(fpu, ['t_', num2str(tstep)]);
    if ~exist(fpu,'dir')
      mkdir(fpu);
    end

    plot_element_data(model, model_data.grid, rb_sim_data.U(:,tstep), plot_params);

    set(gcf,'Color','none');
    set(gca,'Color','none');
    set(gcf,'InvertHardCopy','off');
    showaxes('hide');
    ranges{1} = get(gca,'XLim');
    ranges{2} = get(gca,'YLim');

    export_fig(fullfile(fpu, fn), '-pdf', '-eps', '-png');
    im = export_fig;
    tikzfile(newmodel, fpu, fn, ranges, size(im), 5);
    close(gcf);
  end

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

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

end

%| \docupdate