scripts/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,5]

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_type    = 'nonlinear_richards_trapezoidal';
params.model_type = 'linear_heat_polynomial';
%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
   descr            = richards_fv_descr(params);
   bg_descr         = richards_fv_bg_descr(params);
   [dmodel, rmodel] = gen_models(descr, bg_descr);
   params.mu_ranges = dmodel.mu_ranges;
   params.mu_names  = dmodel.mu_names;
   mu_default = cellfun(@mean, dmodel.mu_ranges);

   model_data = gen_model_data(dmodel);
 case 1 % single detailed simulation and plot
   step1_detailed_simulation;
 case 5 % reduced basis
   detailed_data = gen_detailed_data(rmodel, model_data);
   save(fullfile(rbmatlabresult, detailedfname), 'detailed_data', 'rmodel');
 case 6
   load(fullfile(rbmatlabresult,detailedfname));
   step6_demo_rb_gui;
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