test_subgrid_operators.m

function ok = test_subgrid_operators


  bg_descr.rb_problem_type = 'LinEvol';
  bg_descr.detailed_data_constructor = @Test.DetailedData;
  bg_descr.reduced_data_constructor  = @LinEvol.ReducedData;

  ok = true;

  mm = minimal_model;
  mm.newton_steps = 10;
  mm.newton_epsilon = 1e-10;
  mm.two_phase = TwoPhaseData.Michel;
  mm = mm.two_phase.default_descr(mm);
  mm.mean_ptr        = @(X) harmmean(max(X, 0), 2);
  mm.mean_deriv1_ptr = @(X) 2./( (1./X(:,1) + 1./(X(:,2))).^2 .* X(:,1).^2 );
  mm.mean_deriv2_ptr = @(X) 2./( (1./X(:,1) + 1./(X(:,2))).^2 .* X(:,2).^2 );

  dmodel = TwoPhaseFlow.DetailedModel(mm);
  model_data = gen_model_data(dmodel);
  subfields(1).name = 'saturation';
  subfields(2).name = 'velocity';
  subfields(3).name = 'pressure';
  subfields(1).shortname = 'S';
  subfields(2).shortname = 'U';
  subfields(3).shortname = 'P';
  subfields(1).length = model_data.grid.nelements;
  subfields(2).length = model_data.gn_edges;
  subfields(3).length = model_data.grid.nelements;
  subfields(1).ranges = [0 1];
  subfields(2).ranges = [-2 2];
  subfields(3).ranges = [-5 5];
  mm.subfields = subfields;
  dmodel = TwoPhaseFlow.DetailedModel(mm);
  ok = optest(dmodel, bg_descr, Fv.TwoPhase.SaturationSpace) & ok;
  ok = optest(dmodel, bg_descr, Fv.TwoPhase.VelocitySpace) & ok;
end

function mm = minimal_model
  % declare minimal model
  mm                                  = [];
  mm.name                             = 'minimal_model';
  mm.mu_names                         = {'par1', 'par2', 'par3'};
  mm.mu_ranges                        = {[0 1], [0 2], [0 3]};
  mm.par1                             = 0;
  mm.par2                             = 0;
  mm.par3                             = 0;
  mm.debug                            = 0;
  mm.force_delete                     = 1;
  mm.filecache_ignore_fields_in_model = {'subfields'};
  mm.detailed_simulation              = @(model, model_data) struct('u', model.get_mu(model));
  mm.mass_matrix                      = @fv_mass_matrix;
  mm.get_inner_product_matrix         = @(detailed_data) detailed_data.W;
  mm.l2_error_sequence_algorithm      = @fv_l2_error;

  mm.gridtype                              = 'rectgrid';
  mm.xrange                                = [0 1];
  mm.yrange                                = [0 1];
  mm.xnumintervals                         = 5;
  mm.ynumintervals                         = 5;
  mm.bnd_rect_corner1                      = [-1 -1];
  mm.bnd_rect_corner2                      = [2 2];
  mm.bnd_rect_index                        = -2;

  mm.verbose                               = 9;
  mm.decomp_mode                           = 0;

  mm.new = 1;
  mm.operators_ptr = @(x) x;
  %  mm.init_values_algorithm = @(x) x;

  mm.error_norm = 'l2';
end

function ok = optest(dmodel, bg_descr, ophandle)
  disp(['Testing for operator: ', class(ophandle)]);
  mm = dmodel.descr;

  M_train      = ParameterSampling.Random(30);
  rb_generator = SnapshotsGenerator.Random(dmodel, 'random', [], true, false);
  ei_gen1      = SnapshotsGenerator.SpaceOpEvals(dmodel, 'implicit', rb_generator, ophandle, true, false);

  ei_extension1 = Greedy.Plugin.EI(ei_gen1);
  ei_extension1.stop_Mmax = 200;

  ei_greedy1   = Greedy.Algorithm(ei_extension1, M_train);
  ei_greedy1.stop_timeout = 60*60;
  ei_greedy1.stop_epsilon = 1e-9;

  bg_descr.bg_algorithm = ei_greedy1;
  red_model1     = Test.ReducedModel(dmodel, bg_descr);

  model_data     = gen_model_data(dmodel);
  ok = 1;
  try
    evalc('detailed_data1 = gen_detailed_data(red_model1, model_data);');
  catch exception
    ok = 0;
    disp('ei generation failed');
    disp(getReport(exception));
  end

  TM = detailed_data1.datatree.TM;
  BM = detailed_data1.datatree.BM;
  QM = detailed_data1.datatree.QM;

  [nmd, eind, eind_local] = ophandle.compute_TM_global(model_data, TM);
  reduced_data.grid_local_ext    = nmd.grid_local_ext;
  reduced_data.grid              = nmd.grid_local_ext;
  reduced_data.TM_global_args    = nmd.TM_global_args;
  reduced_data.gEI               = nmd.gEI;
  reduced_data.gn_edges          = nmd.gn_edges;
  reduced_data.gn_inner_edges    = nmd.gn_inner_edges;
  reduced_data.gn_boundary_edges = nmd.gn_boundary_edges;
  reduced_data.TM_global         = eind;
  reduced_data.TM_local          = eind_local;
  reduced_data.BM                = BM;
  reduced_data.opdata = ophandle.fill_opdata(red_model1, detailed_data1.datatree);

  M_train.init_sample(dmodel);
  dmodel.set_mu(M_train.sample(1,:));
  sim_data = rb_generator.generate(dmodel, model_data);

  if isfield(mm, 'subfields')
    [names, ii, jj] = intersect({mm.subfields(:).shortname}, ophandle.arg_vars_short);
    TMg = reduced_data.TM_global_args(jj);
    for i = 1:length(names)
      ei_sim_data.(names{i}) = sim_data.(names{i})(TMg{i});
    end
    offset = size(sim_data.(names{1}), 1);
    for i = 2:length(names)
      TMg{i} = TMg{i} + offset;
      offset = offset + size(sim_data.(names{i}), 1);
    end
    TMg = [TMg{:}];
  else
    ei_sim_data = sim_data(reduced_data.TM_global);
  end

  max_err_sequence = get_field_on_active_child(detailed_data1.datatree, ...
      'max_err_sequence', red_model1);

  ok = ok & (max_err_sequence(end) < 1e-10);
  ok = ok & length(max_err_sequence) <= 31;
  if ok
    disp('epsilon was reached');
  else
    disp('ERROR: epsilon was not reached');
  end

  [full_res,  dummy, full_J]     = ...
    ophandle.apply(mm, model_data, sim_data);
  [subgrid_res dummy, subgrid_J] = ...
    ophandle.apply(mm, reduced_data, ei_sim_data, reduced_data.TM_local);

  oknow = all(abs(full_res(TM) - subgrid_res) < 1e-9);
  if ~oknow
    disp('ERROR: subgrid extraction failed');
    keyboard;
  end
  ok = ok & oknow;
  oknow = all(abs(full_J(TM, TMg) - subgrid_J) < 1e-9);
  if ~oknow
    disp('ERROR: subgrid extraction for jacobian failed');
  end
  sigma = BM \ subgrid_res;
  oknow = all(abs(full_res - QM * sigma) < 1e-9);
  if ~oknow
    disp('ERROR: EI reconstruction failed');
    keyboard;
  end
  ok = ok & oknow;
end