twophase_flow.m

function outfile = twophase_flow(steps, params)
  if nargin == 0
    steps = [0, 1];
  end
  if nargin <= 1
    params = [];
  end

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

  if ~isfield(params, 'size')
    params.size = 0.2;
  end
  if ~isfield(params, 'mode')
    params.mode = 'normal';
  end

  mkdir(fullfile(rbmatlabresult, 'TwoPhaseFlow'));
  outfile_prefix = fullfile(rbmatlabresult, 'TwoPhaseFlow', params.id, 'Script');
  imsavepath = fullfile(getenv('BASEDIR'), 'images', 'TwoPhaseFlow');
  [dummy, sha] = system('git log --no-color --format=''%H'' HEAD | head -1');
  outfile_elems = {'params', 'sha'};

  for si = 1:length(steps)

    step = steps(si);

    switch(step)
      case 0
        descr = twophase_descr(params);

        dmodel = gen_detailed_model(descr);

        model_data = gen_model_data(dmodel);

        bg_descr.rb_problem_type = 'TwoPhaseFlow';
        bg_descr.train_sample_mode = 'single';
        bg_descr.pca = true;

        rmodel = gen_reduced_model(dmodel, bg_descr);

        outfile_elems = [ outfile_elems, {'descr', 'dmodel', 'model_data', 'rmodel'}];
        outfile = [outfile_prefix, 'Step0'];
        save(outfile, outfile_elems{:});
      case 1
        oldoutfile = [outfile_prefix, 'Step0'];
        load(oldoutfile);
        outfile_elems = unique([outfile_elems, matfile_fields(oldoutfile)]);

        tic;
        sim_data = detailed_simulation(dmodel, model_data);
        t = toc;
        sim_data.tictoc = t;

        outfile_elems = [ outfile_elems, {'sim_data'}];
        outfile = [outfile_prefix, 'Step1'];
        save(outfile, outfile_elems{:});

      case 2
        if exist([outfile_prefix, 'Step1.mat'], 'file')
          oldoutfile = [outfile_prefix, 'Step1'];
        else
          oldoutfile = [outfile_prefix, 'Step0'];
        end
        load(oldoutfile);
        outfile_elems = unique([outfile_elems, matfile_fields(oldoutfile)]);
        if exist([outfile_prefix, 'Step1'], 'file')
          load([outfile_prefix, 'Step1']);
        end

        detailed_data = gen_detailed_data(rmodel, model_data);

        outfile_elems = [ outfile_elems, {'detailed_data'}];
        outfile = [outfile_prefix, 'Step2'];
        save(outfile, outfile_elems{:});
      case 3
        oldoutfile = [outfile_prefix, 'Step2'];
        load(oldoutfile);
        outfile_elems = unique([outfile_elems, matfile_fields(oldoutfile)]);
        reduced_data = gen_reduced_data(rmodel, detailed_data);
        tic;
        rb_sim_data = rb_simulation(rmodel, reduced_data);
        t = toc;
        rb_sim_data.tictoc = t;

        tic;
        rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);
        t = toc;
        rb_sim_data.rtictoc = t;

        plot_sim_data(dmodel, model_data, rb_sim_data, []);

        if exist('sim_data', 'var')
          l2_errs     = dmodel.l2_error_sequence_complete(rb_sim_data, sim_data, detailed_data.model_data);
          linfty_errs = dmodel.linfty_error_sequence_complete(rb_sim_data, sim_data, detailed_data.model_data);

          figure;
          subplot(1,2,1); plot_error_sequence(dmodel, l2_errs);     title('L^2-error');
          subplot(1,2,2); plot_error_sequence(dmodel, linfty_errs); title('L^{\infty}-error');

          data = [0:rmodel.descr.dt:rmodel.descr.T; ...
            l2_errs.S; linfty_errs.S; ...
            l2_errs.U; linfty_errs.U; ...
            l2_errs.P; linfty_errs.P];
          datafile = fullfile(imsavepath, 'errors_over_time.dat');

          print_datatable(datafile, {'t', 'sat_l2', 'sat_linfty', 'vel_l2', 'vel_linfty', 'prs_l2', 'prs_linfty'}, data);

          tmpcell = cellfun(@(X) rb_sim_data.(X) - sim_data.(X), {'S', 'U', 'P'}, 'UniformOutput', false);
          diff_data_rb=cell2struct(tmpcell, {'S', 'U', 'P'}, 2);
          plot_sim_data(dmodel, model_data, diff_data_rb, []);

          datafile = fullfile(imsavepath, 'basis_gen.dat');
          N = reduced_data.N;
          Ns = get_by_description(N, 'saturation');
          Nu = get_by_description(N, 'velocity');
          Np = get_by_description(N, 'pressure');
          M = reduced_data.M;
          Ms = get_by_description(M, 'L_saturation');
          Mu = get_by_description(M, 'L_velocity');
          l2_errs.tot = l2_errs.S + l2_errs.U + l2_errs.P;
          linfty_errs.tot = linfty_errs.S + linfty_errs.U + linfty_errs.P;
          data = [ Ns, Nu, Np, Ms, Mu, ...
            model_data.grid.nelements, ...
            max(l2_errs.S), max(linfty_errs.S), ...
            max(l2_errs.U), max(linfty_errs.U), ...
            max(l2_errs.P), max(linfty_errs.P), ...
            max(l2_errs.tot), max(linfty_errs.tot), ...
            sim_data.tictoc, rb_sim_data.rtictoc]';

          print_datatable(datafile, {'Ns', 'Nu', 'Np', 'Ms', 'Mu', 'H', 'sat_l2', 'sat_linfty', 'vel_l2', 'vel_linfty', 'prs_l2', 'prs_linfty', 'l2', 'linfty', 'time', 'rtime'}, data);
        end

        outfile_elems = [ outfile_elems, {'reduced_data', 'rb_sim_data'}];
        outfile = [outfile_prefix, 'Step3'];
        save(outfile, outfile_elems{:});

    end
  end

end

function fnames = matfile_fields(matfile)
  temp = whos('-file', matfile);
  fnames = arrayfun(@(X) X.name, temp, 'UniformOutput', false)';
end