laplace_dune_descr.m

Go to the documentation of this file.

function descr = laplace_dune_descr(params)
% function model = convdiff_dune_model(params)
% convection diffusion example using dune-rb for high dimensional computations
% and data structures

descr = LinStat.descr_default;

descr.name            = 'laplace_dune';
descr.rb_problem_type = 'LinStatDune';
descr.RB_detailed_data_constructor = @LinStat.DetailedData;
%descr.RB_reduced_data_constructor  = @LinEvol.ReducedData;
%% data that is definitely used outside of the detailed simulations
%% These values are now set in dunerbconvdiff( 'init_model')

if (~isempty(getenv('DUNERBHOME')))
%  addpath( fullfile( getenv('DUNERBHOME'), 'mexclient' ) );
  addpath( getenv('DUNERBHOME') );
end

descr.mexptr = @mexclient;

[~] = mexclient('echo', 1);

is_connected = descr.mexptr('is_connected');

if (~is_connected)
  error('No connection to server!');
end

% server_params.serverhost = 'rbevol.uni-muenster.de'
% server_params.port       = '1909'
% mexclient('init_server', server_params)

% copy model data from detailed dune-rb model
tmp_struct = descr.mexptr('init_model');

%model.diffusion = 0;
%model.velocity

descr.mu_names  = tmp_struct.mu_names;
descr.mu_ranges = tmp_struct.mu_ranges;
descr.mu        = descr.mexptr('get_mu');
descr.run_dir   = tmp_struct.cwd;

tmp = descr.mexptr('rb_symbolic_coefficients', {'Laplace'}, {'RHS'});
fns = fieldnames(tmp);
sym_coeffs = [];
for i = 1:length(fns);
  sym_coeffs.(fns{i}) = '@(mu) [ ';
  for j = 1:size(tmp.(fns{i}),2);
    if j > 1
      sym_coeffs.(fns{i}) = [ sym_coeffs.(fns{i}), ', ' ];
    end
    sym_coeffs.(fns{i}) = [ sym_coeffs.(fns{i}), tmp.(fns{i}){j} ];
  end
  sym_coeffs.(fns{i}) = [ sym_coeffs.(fns{i}) , ' ]' ];
end

coeff_ops = structfun(@eval, sym_coeffs, 'UniformOutput', false);

descr.coeff_ops.LI = coeff_ops.Laplace;
descr.coeff_ops.f  = coeff_ops.RHS;
descr.L_I_inv_norm_bound = @(mu) 1/mu(1);

descr.filecache_ignore_fields_in_model = {'coeff_ops', 'mu'};

descr.verbose = 1;
descr.mu_list = {[5 0.5], [1 0.5], [1 0]};
descr.debug = 0;

descr.data_const_in_time     = 1;

descr.RB_train_num_intervals = 5;

%descr.RB_train_size          = [10 10];
descr.RB_train_size          = 200;

descr.RB_stop_Nmax           = 200;

descr.RB_val_size            = 10;

descr.RB_stop_epsilon        = 3e-5;

descr.RB_refinement_mode     = 'adaptive';

descr.RB_customized_basis_generation_ptr = @my_greedy_basisgen;
%% choose generation mode
%descr.RB_generation_mode = 'greedy_uniform_fixed';
%% alternative:
%% RB_generation_mode = 'random_fixed';

%% choose extension method
%descr.RB_extension_algorithm = @dune_RB_extension_PCA_fixspace;
%% stop timeout
%descr.RB_stop_timeout = 60*60; % 1 hour
%% stop on maximum number of base functions
%descr.RB_stop_Nmax = 30;
%% stop on epsilon
%descr.RB_stop_epsilon = 1e-9; % 0.003 is already realized by init data!!
%% stop on overfitting (needs a validation test set)
%descr.RB_stop_max_val_train_ratio = inf;

%%descr.RB_test_size = 1000;

%% fix the rand seed in case of RB_generation_mode = random_fixed
%descr.RB_train_rand_seed = 0.5;
%% trainig set size in case of RB_generation_mode = random_fixed
%descr.RB_train_size = 50;

%% choose the error indicator
%descr.RB_error_indicator = 'estimator'; % Delta from rb_simulation
%% alternative:
% model.RB_error_indicator = 'error'; % true error

function RB = my_basisgen(rmodel, sdd)

  dmodel = rmodel.detailed_model;
  descr = dmodel.descr;
  if isempty(descr.mu_list)
    descr.mu_list = {cellfun(@mean, rmodel.descr.mu_ranges, 'UniformOutput', true)};
  end

  sdd.mu_list = descr.mu_list;
  RB = [];

  for m = 1:length(descr.mu_list)
    next = descr.mexptr('rb_extension_PCA', descr.mu_list{m}, 1);
    if next == 0
      continue;
    end
    RB = [ RB, next ];
  end
  
function RB = my_greedy_basisgen(rmodel, sdd)
disp('entered my_Greedy_basisgen');

dmodel = rmodel.detailed_model;
set_mu(dmodel, [1 0]);
set_mu(rmodel, [1 0]);
descr = rmodel.detailed_model.descr;
sdd.RB = descr.mexptr('init_data_basis');

reduced_data = gen_reduced_data(rmodel, sdd);

% Kommentar MD: Alle Description Felder, die mit RB_ beginnen, werden von
% gen_reduced_model in das Feld bg_descr kopiert (ohne den Präfix 'RB_')
% (bg_descr.x == descr.RB_x)
bg_descr = rmodel.bg_descr;

%ps = ParameterSampling.Uniform(bg_descr.train_size);
%ps.init_sample(dmodel);
%mus = ps.sample';
mus = rand_uniform(bg_descr.train_size,rmodel.mu_ranges);
disp(mus);
max_err_est = 1e10;
mu_max_list = [];
mu_list = [];
while( max_err_est> bg_descr.stop_epsilon) && ...
      (size(sdd.RB,2)< bg_descr.stop_Nmax)
  max_err_est = 0;
  mu_max = [];
  errs = zeros(1, size(mus,2));
  for i = 1:size(mus,2);
    set_mu(rmodel, mus(:,i));
    rb_sim_data = rb_simulation(rmodel,reduced_data);
    errs(i) = rb_sim_data.Delta;
    if rb_sim_data.Delta > max_err_est
      mu_max = mus(:,i);
      max_err_est = rb_sim_data.Delta;
    end
    
  end
  disp(errs);
  mu_max_list = [mu_max_list, max_err_est];
  mu_list = [mu_list, mu_max];
  disp(['max error estimator: ',num2str(max_err_est),...
        ' for mu = ',num2str(mu_max')]);
  if (isempty(mu_max))
    break;
  end
  next = descr.mexptr('rb_extension_PCA', mu_max, 1);
  sdd.RB = [sdd.RB, next];

  reduced_data = gen_reduced_data(rmodel,sdd);
  disp(['new basis size: ',num2str(size(sdd.RB,2))]);
end

save('/tmp/maxes', 'mu_list', 'mu_max_list');

RB = sdd.RB;