1 function [descr, rdescr, dmodel, rmodel] =
newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
2 %
function [descr, rdescr, dmodel, rmodel] =
newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
3 % small script demonstrating a buckley leverett problem
4 % discretization and RB model reduction
6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 %%%%%% Select here, what is to be performed with the richards data %%%%%%%%
8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 %step = 1 % single detailed simulation with given data and plot. Run
10 %
this with varying parameters mu until sure that scheme
11 % is stable. Modify dt or the data-functions accordingly,
12 % until a nice parameter-domain with uniformly stable
13 % detailed scheme is obtained.
14 %step = 2 % generate colateral reduced basis of L_E
operator
15 %step = 3 % compare single detailed simulation with and without
16 % interpolated space
operator
17 %step = 4 % generate dummy reduced basis from single trajectory and check,
if
18 % ei_interpolation with projection on
this space maintains
19 % result. A simple reduced simulation can also be
20 % performed. All results should be visually identical
21 %step = 5 % generate reduced basis
22 %step = 6 % time measurement of reduced simulation and
23 % use reduced basis in rb_demo_gui
25 % written by M. Drohmann
27 %steps = [0,2,3,4,5,6]
35 noof_ei_extensions = 1;
37 model_size =
'normal';
46 warning(
'off',
'RBMatlab:rb_simulation:runtime');
47 warning(
'off',
'MATLAB:cholinc:ArgInfToBeRemoved');
49 imsavepath = fullfile(rbmatlabresult,
'images');
53 for si=1:length(steps)
57 params.use_laplacian = use_laplacian;
59 params.model_type = 'exponential_diffusion';
60 extrapar.Mstrich = 20;
62 params.model_size = model_size;
63 %params.model_type = 'eoc_nonlinear';
64 %params.model_size = 'big';
65 %params.model_size = 'medium';
66 % params.model_size = 'small';
67 % params.model_type = 'nonlinear';
69 params.step7_outputfile = fullfile(rbmatlabresult, ['newton_', params.model_type, '_step7.mat']);
70 params.step8_outputfile = fullfile(rbmatlabresult, ['newton_', params.model_type, '_step8.mat']);
77 % email to which status messages about processing mu vectors is sent
78 % params.email = 'mdrohmann@uni-muenster.de';
80 %% parameters for visualization
81 plot_params.show_colorbar = 1;
82 plot_params.colorbar_mode = 'EastOutside';
83 plot_params.no_lines = true;
84 plot_params.plot_type = 'contour';
85 % plot_params.clim = [0, 0.5];
86 % plot_params.clim_tight = true;
88 plot_params.bind_to_model = true;
89 plot_params.yscale_uicontrols = 0.6;
91 % output-filenames in rbmatlabresult
92 detailedfname = ['newton_', params.model_type, '_', num2str(Mstrich), '_detailed_interpol.mat'];
95 case 0 % initialize model data
96 if ~isempty(getenv('MACHINEFILE'))
97 machinefile = getenv('MACHINEFILE');
98 FakeMPI.MPI_Init(machinefile, fullfile(rbmatlabhome, 'fakeMPI', 'run_parCompErrs.sh'));
102 descr = newton_oo_model(params);
103 rdescr.rb_problem_type = descr.rb_problem_type;
105 rdescr.train_sample_mode = 'random';
106 rdescr.train_num_intervals = 800;
107 rdescr.stop_max_val_train_ratio = inf;
110 rdescr.stop_timeout = Inf;
111 rdescr.stop_Nmax = 300;
112 rdescr.stop_Mmax = 600;
113 rdescr.stop_epsilon = 1e-4;
114 rdescr.val_size = 20;
115 rdescr.stop_max_val_train_ratio = 1.10;
116 rdescr.ei_minimum_residual = 1e-10;
117 rdescr.extension_M_by_N_r = M_by_N_ratio;
118 rdescr.noof_ei_extensions = noof_ei_extensions;
119 rdescr.maximum_temporary_error_growth_factor = 1e3;
120 rdescr.train_num_intervals = [4,1,2];
121 rdescr.max_refinement_level = 8;
122 rdescr.Mmax_small = 30;
123 rdescr.ei_epsilon_small = 1e-3;
126 rdescr.stop_Nmax = 111;
127 rdescr.Mmax_small = rdescr.Mmax_small + Mstrich;
130 rdescr.train_num_intervals = [7,2,4];
131 rdescr.Mmax_small = 0;
133 rdescr.indicator_mode = 'error';
134 rdescr.stop_Nmax = 100;
135 rdescr.stop_Mmax = 426;
136 rdescr.stop_epsilon = 0;
139 rdescr.stop_Nmax = 100;
144 rdescr.refinement_theta = 0.10;
148 rmodel.Mstrich = Mstrich;
150 rmodel.enable_error_estimator = false;
153 dmodel.num_cpus = num_cpus;
154 rmodel.num_cpus = num_cpus;
156 params.mu_ranges = rmodel.mu_ranges;
157 params.mu_names = rmodel.mu_names;
159 model_data = gen_model_data(dmodel);
160 case 1 % single detailed simulation and plot
163 tmppar.model_type = 'eoc';
165 model = newton_model(tmppar);
167 model.newton_solver = 0;
169 model.dt = model.T/model.nt;
171 err = zeros(model.nt+1,7);
174 model_data = gen_model_data(model);
176 sim_data = detailed_simulation(model, model_data);
178 ffe = @exact_function_heat_equation;
180 exact_data.U = zeros(model_data.grid.nelements, model.nt+1);
182 for tt = 1:model.nt+1
183 model.t = (tt-1) * model.dt;
184 exact_data.U(:,tt) = ffe([model_data.grid.CX(:), model_data.grid.CY(:)], model);
188 err(:,i) =
fv_l2_error(sim_data.U, exact_data.U, model_data.grid, model);
190 eoc(i-1) = log(max(err(:,i-1))/max(err(:,i)))/log(2);
191 disp(['Error: ', num2str(err(end,i)), 'EOC: ', num2str(eoc(i-1))]);
195 model.xnumintervals = round(model.xnumintervals.*2);
196 model.ynumintervals = round(model.ynumintervals.*2);
199 tmppar.model_type = 'eoc_nonlinear';
200 tmppar.local_mass_matrix = @fv_local_mass_matrix_rect;
202 tmppar.evaluate_basis = @fv_evaluate_basis;
205 model = newton_model(tmppar);
207 % model.newton_regularisation = 0;
209 model.dt = model.T/model.nt;
210 model.ynumintervals = 1;
211 err = zeros(model.nt+1,5);
214 disp(['Computing EOC step no. ', num2str(i)]);
215 model_data = gen_model_data(model);
217 sim_data = detailed_simulation(model, model_data);
221 width = model_data.grid.X(2) - model_data.grid.X(1);
223 Fe = @(einds, loc, grid, params) ffe([grid.X(einds) + loc(:,1)*width,...
224 grid.Y(einds) + loc(:,2)*width], ...
226 tmppar.diff_k0 = model.diff_k0;
227 exact_data.U = zeros(model_data.grid.nelements, model.nt+1);
228 for tt = 1:model.nt+1
229 model.t = (tt-1) * model.dt;
230 exact_data.U(:,tt) = ffe([model_data.grid.CX(:), ...
231 model_data.grid.CY(:)], model);
232 % exact_data.U(:,tt) = l2project(Fe, 2, model_data.grid, tmppar);
236 err(:,i) =
fv_l2_error(sim_data.U, exact_data.U, model_data.grid, model);
238 eoc(i-1) = log(max(err(:,i-1))/max(err(:,i)))/log(2);
239 disp(['Error: ', num2str(err(end,i)), ' EOC: ', num2str(eoc(i-1))]);
242 model.xnumintervals = model.xnumintervals.*2;
244 model.dt = model.T / model.nt;
245 %model.ynumintervals = model.ynumintervals.*2;
248 case 2 % empirical interpolation of space operator
250 case 3 % compare detailed simulation with and without interpolation
251 model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
253 case 4 % construct dummy reduced basis by single trajectory and simulate
254 model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
256 case 5 % reduced basis
257 % load(fullfile(rbmatlabresult,CRBfname));
259 detailed_data = gen_detailed_data(rmodel, model_data);
260 save(fullfile(rbmatlabresult, detailedfname), 'detailed_data', 'rmodel');
262 tmp=load(fullfile(rbmatlabresult,detailedfname));
263 detailed_data=tmp.detailed_data;
266 case 7 % training-error landscape & time comparisons
267 load(fullfile(rbmatlabresult,detailedfname));
269 basis_gen.Mstrich = 0;
270 % number of reduced basis sizes to be tested.
272 % number of collateral reduced basis sizes to be tested.
274 % sample of values for the coupling variable 'c' for which errors are
279 % axis description for plots with coupling variable 'c'
280 cdescr = ['(N,M) = c * (Nmax, ',num2str(M_by_N_ratio), '*Mmax)'];
281 % size of mu vector test set.
284 error_plots = {
'train_set_coupled', ...
291 load(params.step7_outputfile);
292 for i = 1:length(output)
296 disp(
'warning: takes a few hours!');
297 load(fullfile(rbmatlabresult,detailedfname));
301 % number of reduced basis sizes to be tested.
303 % number of collateral reduced basis sizes to be tested.
305 % Mstrich values to be tested
306 Mstrich_samples = [1,2,3,5,8,10,15,20];
307 % sample of values
for the coupling variable
'c' for which
308 % estimators are computed.
310 % size of mu vector test set
312 % axis description
for plots with coupling variable
'c'
313 cdescr =
'(N,M) = c * (Nmax, Mmax)';
317 load(params.step8_outputfile);
318 for i = 1:length(output)
319 output{i}.bound = 1500;
320 output{i}.stab_limit = 10000;
324 load(fullfile(rbmatlabresult,detailedfname));
325 model_data = detailed_data;
326 testdir =
'richards_test_10';
328 model.Nmax = size(detailed_data.RB,2);
330 % generate test data
if not existent
331 rand(
'state',123456);
332 M = rand_uniform(10,params.mu_ranges);
335 % load demo_nonlin_evol_detailed_data_LE_on_RB;
336 reduced_data = gen_reduced_data(model,detailed_data);
338 settings = load(fullfile(rbmatlabtemp,...
339 testdir,
'settings.mat'));
343 % storage
for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
344 errs = zeros(Msamples,Nsamples);
345 model.RB_error_indicator='ei_estimator_test';
347 Mstrich_set = 1:10:100;%[1,2,3,4,5,10,15,20];
349 deltas = zeros(size(M,2),length(Mstrich_set));
350 for mu_ind = 1:size(M,2);
351 disp(['processing parameter vector ',num2str(mu_ind),...
352 '/',num2str(size(M,2))]);
353 for msi = 1:length(Mstrich_set)
355 model = model.set_mu(model,settings.M(:,mu_ind));
361 mtemp.M = extrapar.M;
362 mtemp.Mstrich = Mstrich_set(msi);
363 mtemp.N = extrapar.N;
364 simulation_data = rb_simulation(mtemp,reduced_data);
367 deltas(mu_ind,msi) = simulation_data.Delta(end);
368 simulation_data = rb_reconstruction(mtemp,detailed_data,simulation_data);
369 error = model.error_algorithm(simulation_data.U, sim_data.U, ...
370 detailed_data.grid, mtemp);
371 disp(['true error: ', num2str(max(error))]);
373 %
plot_sim_data(model,detailed_data,simulation_data,plot_params);
376 plot(Mstrich_set,deltas(mu_ind,:));
381 % tmp = load(fullfile(rbmatlabresult, 'exponential_laplace', 'newton_exponential_diffusion_detailed_interpol_new'));
382 tmp = load(fullfile(rbmatlabresult, detailedfname));
383 detailed_data = tmp.detailed_data;
384 mu_set = {[1,0.01,0.2],[1,0.01,0.0],[4,0.01,0.2],[4,0.01,0.0]};
386 plot_params.plot_type =
'contour';
387 plot_params.contour_lines = 7;
388 plot_params.line_width = 2;
389 timeinstants = [0.1,1.0]*model.T;
390 boxed_snapshots =
true;
393 %tmp = load(fullfile(rbmatlabresult,
'strange.mat'));
395 tmp = load(fullfile(rbmatlabresult,detailedfname));
396 detailed_data = tmp.detailed_data;
398 plot_params.no_lines = 1;
400 plot_params.no_axes = 1;
401 plot_params.transparent_background = 1;
402 plot_params.show_colorbar = 0;
403 %plot_params.shrink_factor = 1.13;
404 %plot_params.colormap = Cmap;
405 plot_params.clim = [0.1,0.6];
407 reduced_data = gen_reduced_data(model, detailed_data);
408 model.N = size(detailed_data.RB,2);
409 model.M = size(detailed_data.QM,2);
410 reduced_data = extract_reduced_data_subset(model, reduced_data);
412 cparams.range = cell(length(model.mu_ranges)+1,1);
413 cparams.range{1} = [1, model.nt];
414 cparams.range(2:end) = model.mu_ranges(:);
415 cparams.numintervals = 2*ones(length(model.mu_ranges)+1,1);
418 parammat = cgrid.vertex;
419 parammat(:,1) = round(parammat(:,1));
422 for i=1:length(parammat)
423 tstep = parammat(i,1);
424 newmodel = model.set_mu(model, parammat(i,2:end));
425 rb_sim_data = rb_simulation(newmodel, reduced_data);
426 rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
428 % transformed version
429 plot_element_data(model_data.grid, rb_sim_data.U(:,tstep), plot_params);
430 Uline = rb_sim_data.U(subline,tstep);
431 mustring = strrep(strrep(strrep(strrep(mat2str(parammat(i,2:end)),' ', '_'),'[',''),']',''),'.','p');
433 fp = fullfile(imsavepath, params.model_type);
437 fp = fullfile(fp, ['t_', num2str(tstep)]);
441 fn = ['sample_mu_', mustring];
442 disp(['Processing ', fn, ' ...']);
444 tikzparams.filename = fn;
445 tikzparams.filepath = fp;
446 tikzparams.save_colorbar = 1;
454 print_datatable(fullfile(fp, datafn), subline_title, subline_values);
457 load('richards_affine_detailed_interpol');
459 rdata = gen_reduced_data(model, detailed_data);
461 model.N = size(detailed_data,2);
464 rdata2 = extract_reduced_data_subset(model, rdata);
466 rdata1 = extract_reduced_data_subset(model, rdata);
468 TM1 = rdata1.TM_local{1}(1:10);
469 TM2 = rdata2.TM_local{1}(1:10);
471 gl1 = rdata1.grid_local_ext{1};
472 gl2 = rdata2.grid_local_ext{1};
474 if(any(gl1.CX(TM1) ~= gl2.CX(TM2)))
477 if(any(gl1.CY(TM1) ~= gl2.CY(TM2)))
480 if(any(any(gl1.X(gl1.VI(TM1,:)) ~= gl2.X(gl2.VI(TM2,:)))))
483 if(any(any(gl1.Y(gl1.VI(TM1,:)) ~= gl2.Y(gl2.VI(TM2,:)))))
486 if(any(any(gl1.SX(TM1,:) ~= gl2.SX(TM2,:))))
489 if(any(any(gl1.SY(TM1,:) ~= gl2.SY(TM2,:))))
492 if(any(gl1.EL(TM1,:) ~= gl2.EL(TM2,:)))
495 if(any(gl1.DC(TM1,:) ~= gl2.DC(TM2,:)))
498 if(any(gl1.NX(TM1,:) ~= gl2.NX(TM2,:)))
501 if(any(gl1.NY(TM1,:) ~= gl2.NY(TM2,:)))
504 if(any(gl1.ECX(TM1,:) ~= gl2.ECX(TM2,:)))
507 if(any(gl1.ECY(TM1,:) ~= gl2.ECY(TM2,:)))
511 NBI1 = gl1.NBI(TM1,:);
512 NBI2 = gl2.NBI(TM2,:);
513 if(any(gl1.CX(NBI1) ~= gl2.CX(NBI2)))
516 if(any(gl1.CY(NBI1) ~= gl2.CY(NBI2)))
519 if(any(any(gl1.X(gl1.VI(NBI1,:)) ~= gl2.X(gl2.VI(NBI2,:)))))
522 if(any(any(gl1.Y(gl1.VI(NBI1,:)) ~= gl2.Y(gl2.VI(NBI2,:)))))
525 if(any(any(gl1.SX(NBI1,:) ~= gl2.SX(NBI2,:))))
528 if(any(any(gl1.SY(NBI1,:) ~= gl2.SY(NBI2,:))))
531 for i = 1:length(TM1)
534 case 12 % plot some reduced basis functions
535 tmp = load(fullfile(rbmatlabresult,detailedfname));
536 detailed_data = tmp.detailed_data;
539 plot_params.no_lines = 1;
541 plot_params.axis_equal = true;
542 plot_params.plot_type = 'patch';
543 plot_params.no_axes = 1;
544 plot_params.transparent_background = 1;
545 plot_params.show_colorbar = 0;
547 subline = find(detailed_data.grid.CY > 0.6 ...
548 & detailed_data.grid.CY < 0.6+detailed_data.grid.hmin);
549 subline_title = {
'X'};
550 subline_values = grid.CX(subline)
';
553 plot_element_data(detailed_data.grid, detailed_data.RB(:,i), plot_params);
554 Uline = detailed_data.RB(subline,i);
555 subline_title = [ subline_title, ['rb
',num2str(i) ] ];
556 subline_values = [ subline_values; reshape(Uline, 1, length(Uline)) ];
558 fp = fullfile(imsavepath, params.model_type);
562 fn = ['rb_function_
', num2str(i)];
563 disp(['Processing
', fn, ' ...
']);
565 tikzparams.filename = fn;
566 tikzparams.filepath = fp;
567 tikzparams.save_colorbar = 1;
568 tikzparams.width = 3.5;
569 tikzparams.scaleaxis = '';
571 plot_as_tikzfile(model, tikzparams);
575 datafn = 'rb_sublines.dat
';
576 print_datatable(fullfile(fp, datafn), subline_title, subline_values);
580 error('step-number is unknown!
');
function clim = plot_as_tikzfile(model, params)
postprocesses a figure and write out an image and a text file that can be included in TeX documents...
function [ sim_data , tictoc ] = load_detailed_simulation(m, savepath, params)
load single trajectory of previously saved results.
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
function step1_detailed_simulation()
script performing a single detailed simulation and plotting it.
function [ dmodel , rmodel ] = gen_models(ModelDescr descr,BasisGenDescr bg_descr)
generates an IDetailedModel and an IReducedModel instance from description structures.
function plot_error_landscape(output)
plots an output structure generated by stochastic_error_estimation()
A hierarchical cubegrid of arbitrary dimension.
function p = fv_plot(gridbase grid, dofs, params)
routine plotting a single fv function of fv_functions.
function step2_empirical_interpolation()
script constructing a collateral reduced basis space for localized space operators.
function step4_dummy_reduced_basis()
script constructing a reduced basis from a single trajectory and performing a detailed simulation whe...
function step7_error_landscape()
script generating error landscapes by computing the true error of reduced simulations vs...
function res = cubequadrature(poldeg, func, varargin)
integration of function func over reference cube == unit cube. by Gaussian quadrature exactly integra...
function save_detailed_simulations(model, model_data, M, savepath)
perform loop over detailed simulations and save results or check consistency with existing saved resu...
function step8_estimator_landscape()
script generating landscapes plot data by computing the error estimator of reduced simulations for gi...
function step6_demo_rb_gui()
script comparing time for a reduced and a detailed simulation and starting the demonstration GUI ...
function buckley_leverett()
small script demonstrating a buckley leverett problem discretization and RB model reduction ...
function l2_error = fv_l2_error(U1, U2, W)
function computing the l2-error between the two fv-functions or function sequences in U1...
function [ descr , rdescr , dmodel , rmodel ] = newton(steps, combined, M_by_N_ratio, noof_ei_extensions, use_laplacian, model_size, random, num_cpus, Mstrich)
small script demonstrating a buckley leverett problem discretization and RB model reduction ...
function res = exact_function_plaplace(glob, params)
This is the first function from http://eqworld.ipmnet.ru/en/solutions/npde/npde1201.pdf.
function step3_detailed_ei_simulation()
script performing a detailed simulation with empirical interpolated operators comparing the result wi...
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
function step10_plot_trajectories()
script generating a tikz graphic showing trajectories for certain selected parameters ...