1 % script demonstrating the burgers equation with
explicit fv
2 % discretization, empirical interpolation and RB model reduction
4 % a step needs to be selected in
this script,
default is 1.
6 % example is meant to demonstrate automatic space dimension reduction
8 % The interpolation basis functions are nicely vertical stripes.
9 % The error convergence in ei however is constant until 100 basis
10 % functions then drops to zero. Interesting fact, but clearly M not
11 % variable afterwards!!!
12 % As a reduced basis, two shock-trajectories are used.
14 % Bernard Haasdonk 16.6.2008
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17 %%%%%% Select here, what is to be performed with the burgers data %%%%%%%%%%
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19 step = 1 % single detailed simulation with given data and plot. Run
20 %
this with varying parameters mu until sure that scheme
21 % is stable. Modify dt or the data-functions accordingly,
22 % until a nice parameter-domain with uniformly stable
23 % detailed scheme is obtained.
24 %step = 2 % generate colateral reduced basis of L_E
operator
25 % takes still 20 minutes
if snapshots are already computed!!
26 %step = 3 % compare single detailed simulation with and without
27 % interpolated space
operator
28 %step = 4 % generate dummy reduced basis from trajectories and check,
if
29 % ei_interpolation with projection on
this space maintains
30 % result. A simple reduced simulation can also be
31 % performed. All results should be visually identical
32 %step = 5 % generate reduced basis
33 %step = 6 % time measurement of reduced simulation
35 %step = 7 % generate error-landscape over varying N (and eventually M)
36 % can take several hours!!!
37 %step = 8 % generate different Figures
for presentation and time measurement
39 %step = 10 % start
demo_rb_gui without error estimator
41 % output-filenames in rbmatlabtemp
42 CRBfname =
'riemann_burgers_CRB.mat';
43 detailedfname =
'riemann_burgers_detailed.mat';
44 testdir =
'riemann_burgers_test_10';
45 results_fname_base =
'riemann_burgers_step7_result';
46 ei_detailed_savepath =
'riemann_burgers_ei_detailed';
47 ei_operator_savepath =
'riemann_burgers_detailed';
49 % old settings, now in riemann_burgers_model:
51 model = riemann_burgers_model([]);
55 %% rectangular or triangular
56 %%params.gridtype =
'triagrid';
57 %%params.grid_initfile =
'2dsquaretriang';
59 %params.gridtype =
'rectgrid';
60 %params.xnumintervals = 100;
61 %params.ynumintervals = 100;
62 %params.xrange = [0,1];
63 %params.yrange = [0,1];
65 % set upper and lower to neuman noflow
66 % set left and right to dirichlet
69 %params.bnd_rect_corner1 = [-1, 0-eps, 1-eps; ...
71 %params.bnd_rect_corner2 = [ 2, 0+eps, 1+eps; ...
74 % -1 means dirichlet, -2 means neuman
75 %params.bnd_rect_index = [-2,-1,-1];
79 %params.nt = 100; % guess and adjust later
82 %params.nt = 80; % guess and adjust later
83 %params.nt = 100; % guess and adjust later
85 %params.nt = 50; % guess and adjust later
86 %params.400 => non-bounded u
90 %params.axis_equal = 1;
91 %params.clim = [-1,1];
94 %params.name_init_values =
'blobs';
97 %params.beta = 1; % init data weight between 0 and 1
99 % name of
function in rbmatlab/datafunc/init_values
100 %params.name_init_values =
'waveproduct';
101 %params.name_init_values =
'leftright';
102 %params.name_init_values =
'as_dirichlet';
104 %params.c_init_max = 1.0;
105 %params.c_init_phase_x = 0.0;
106 %params.c_init_phase_y = 0.0;
107 %params.c_init_freq_x = 2*pi;
108 %params.c_init_freq_y = 2*pi;
110 %params.c_init_lo = 1.0;
111 %params.c_init_lo = 0.0;
112 %params.c_init_lo = -1.0;
115 %params.name_dirichlet_values =
'leftright';
116 %params.c_dir_left = -1.0;
117 %params.c_dir_right = 0.0;
118 %params.dir_middle = 0.5;
120 % name of
function in RBmatlab/datafunc/neuman_values
121 %params.name_neuman_values =
'zero';
122 %params.name_neuman_values =
'homogeneous';
126 %params.name_flux =
'burgers_parabola';
127 %params.name_flux =
'burgers';
128 %
for simplifying computation in
case of linear flux:
129 %params.flux_linear = 0;
130 % horizontal right velocity field
133 %params.flux_quad_degree = 2;
134 %params.flux_pdeg = 2; % burgers exponent
137 %params.mu_names = {
'c_init_lo',
'vrot_angle'};
138 % themiddle parameter makes u_init not decomposable...
139 %params.mu_names = {
'c_dir_left',
'c_dir_right',
'dir_middle'};
140 %params.mu_names = {
'c_dir_left',
'c_dir_right',
'flux_vx'};
141 %params.mu_ranges = {[0,1],[0,1],[0,1]};
142 %params.mu_ranges = {[-1,1],[-1,1],[-1,1]};
143 %params.rb_problem_type =
'nonlin_evol';
144 %params.problem_type =
'lin_evol';
146 % finite volume settings
147 %params.name_diffusive_num_flux =
'none';
148 %params.name_convective_num_flux =
'lax-friedrichs';
149 %params.lxf_lambda = 1.66
150 %params.lxf_lambda = 1.25;
151 %params.lxf_lambda = 0.35355;
152 %params.lxf_lambda = fv_search_max_lxf_lambda([],params);
154 %params.name_convective_num_flux =
'engquist-osher';
156 % projection of analytical initial data to discrete
function
157 %params.init_values_algorithm =
'fv_init_values';
158 %params.implicit_operators_algorithm =
'fv_operators_implicit';
159 %
get mass matrix
for inner product computation from DOF-vectors
160 %params.inner_product_matrix_algorithm =
'fv_inner_product_matrix';
161 %params.data_const_in_time = 1;
163 % settings
for empirical interpolation
164 %params.L_E_local_name =
'fv_conv_explicit_space';
165 %par.range = params.mu_ranges;
166 %params.ei_numintervals = [4,4,4]; % 5^3 parameter grid
167 %params.ei_numintervals = [2,2,2]; % 5^3 parameter grid
168 %params.ei_numintervals = [4,4]; % 5^2 parameter grid
171 %params.ei_detailed_savepath = ei_detailed_savepath;
172 %params.ei_operator_savepath = ei_operator_savepath;
173 %params.ei_time_indices = 1:params.nt;
174 %params.CRB_generation_mode =
'param-time-space-grid';
175 %params.ei_target_error =
'interpol';
177 % settings
for reduced basis generation
178 %params.detailed_simulation_algorithm =
'detailed_simulation';
179 %params.operators_algorithm =
'fv_operators_implicit_explicit';
180 %params.init_values_algorithm =
'fv_init_values';
181 %params.inner_product_matrix_algorithm =
'fv_inner_product_matrix';
182 %params.error_algorithm =
'fv_error';
183 %params.lxf_lambda = 1.0194e+003;
184 %params.data_const_in_time = 1;
185 %params.error_norm =
'l2';
186 %params.RB_extension_algorithm =
'RB_extension_PCA_fixspace';
187 %
'RB_extension_max_error_snapshot'};
188 %params.RB_stop_timeout = 2*60*60; % 2 hours
189 %params.RB_stop_epsilon = 1e-5;
190 %params.RB_error_indicator =
'error'; %
true error
191 %params.RB_error_indicator =
'estimator'; % Delta from rb_simulation
192 %params.RB_stop_Nmax = 100;
193 %params.RB_generation_mode =
'uniform_fixed';
194 %params.RB_generation_mode =
'PCA_trajectories';
195 %params.RB_mu_list = {[1,0,1],[0,1,-1]}; % two shock-trajectories in
196 % different direction
197 %params.RB_numintervals = params.ei_numintervals;
198 %params.RB_detailed_train_savepath = params.ei_detailed_savepath;
201 %params.yscale_uicontrols = 0.6;
202 %params.xscale_gui = 0.5;
204 model_data = gen_model_data(model);
207 case 1 % single detailed simulation and plot
208 disp(
'performing single detailed simulation')
209 % grid = construct_grid(model);
210 % shock moving to left:
211 % model = set_mu(model,[-0.5,1,-1]); % cdir_left, c_dir_right, vx
213 % model = set_mu(model,[-1,1,-1]); % cdir_left, c_dir_right, vx
214 % symmetric rarefaction wave:
215 % model = set_mu(model,[1,-1,-1]); % cdir_left, c_dir_right, vx
216 % moving shock with velocity 0.5:
217 model = set_mu(model,[0,1,-1]); % cdir_left, c_dir_right, vx
218 sim_data = detailed_simulation(model, model_data);
221 case 2 % empirical interpolation of space
operator
223 disp(
'constructing collateral reduced basis')
224 model.RB_generation_mode = 'none';
225 detailed_data = gen_detailed_data(model,model_data);
226 save(fullfile(rbmatlabtemp,CRBfname),...
227 'detailed_data','model');
228 detailed_data.RB = zeros(detailed_data.grid.nelements,1);
230 plot_params.plot = model.plot;
231 plot_detailed_data(model,detailed_data,plot_params);
234 disp(['computation time for collateral basis: ',num2str(t)]);
235 case 3 % compare detailed simulation with and without interpolation
236 disp(['comparison between detailed simulation with and without', ...
238 tmp = load(fullfile(rbmatlabtemp,CRBfname));
239 detailed_data = tmp.detailed_data;
240 % params.flux_pdeg = 2;
241 sim_data = detailed_simulation(model,model_data);
242 plot_params.title = 'detailed simulation without interpolation';
244 model.M = size(detailed_data.QM{1},2);
245 % model.M = size(detailed_data.QM,2);
247 sim_data_interpol = nonlin_evol_detailed_ei_simulation(model,detailed_data);
248 plot_params.title =
'detailed simulation with empirical interpolation';
249 plot_sim_data(model,model_data,sim_data_interpol,plot_params);
250 plot_params.title =
'error';
252 err.U = err.U-sim_data_interpol.U;
254 disp([
'maximum l-infty error:',num2str(max(max(abs(err.U))))]);
256 case 4 % construct dummy reduced basis by single trajectory and simulate
257 load(fullfile(rbmatlabtemp,CRBfname));
258 model = riemann_burgers_model;
259 % set mu values
new, as may be different in stored file
260 model.c_dir_left = 1.0;
261 model.c_dir_right = 0.0;
262 model.dir_middle = 0.5;
264 disp(
'detailed interpolated simulation for basis construction:')
265 model.M = size(detailed_data.QM{1},2);
266 sim_data_interpol = ...
267 nonlin_evol_detailed_ei_simulation(model,detailed_data);
268 % A = feval(params.inner_product_matrix_algorithm, ...
269 % detailed_data.grid,params);
270 UON = model.orthonormalize(model,model_data,sim_data_interpol.U,1e-4);
271 % UON = model.orthonormalize(model,model_data,sim_data_interpol.U,1e-4);
272 detailed_data.RB = UON;
273 model.N = size(detailed_data.RB,2);
274 model.M = size(detailed_data.QM{1},2);
275 save(fullfile(rbmatlabtemp,
'riemann_burgers_temp'),
'model',
'detailed_data')
277 disp(
'reduced simulation:')
278 reduced_data = gen_reduced_data(model,detailed_data);
279 sim_data = rb_simulation(model,reduced_data);
280 sim_data = rb_reconstruction(model,detailed_data,sim_data);
282 disp('detailed interpolated and rb-projected simulation:')
283 model.M = size(detailed_data.QM{1},2);
284 model.N = size(detailed_data.RB,2);
285 sim_data_ei_rb_proj = ...
286 nonlin_evol_detailed_ei_rb_proj_simulation(model,detailed_data);
288 plot_params.title =
'reduced simulation result';
290 plot_params.title =
'detailed ei_interpol simulation result';
291 plot_sim_data(model,model_data,sim_data_interpol,plot_params);
292 plot_params.title =
'detailed ei_interpol and rb_projected simulation result';
293 plot_sim_data(model,model_data,sim_data_ei_rb_proj,plot_params);
295 case 5 % reduced basis construction
296 disp(
'constructing reduced basis')
297 tmp = load(fullfile(rbmatlabtemp,CRBfname));
298 model = riemann_burgers_model;
299 detailed_data = tmp.detailed_data;
301 % set these fields such that
demo_rb_gui works nicely
302 model.N = size(detailed_data.RB,2);
303 model.M = size(detailed_data.QM{1},2);
304 save(fullfile(rbmatlabtemp,detailedfname),...
305 'detailed_data',
'model');
307 if isfield(detailed_data.RB_info,
'max_err_sequence')
308 plot(detailed_data.RB_info.max_err_sequence);
309 set(gca,'Yscale','log');
310 title('RB-generation error convergence');
312 disp(['skipping plot of error decrease as basis generation does' ...
316 case 6 % time measurement of reduced simulation
318 load(fullfile(rbmatlabtemp,detailedfname));
319 disp('reduced simulation:')
320 model.N = size(detailed_data.RB,2);
321 model.M = size(detailed_data.QM{1},2);
322 reduced_data = gen_reduced_data(model,detailed_data);
324 sim_data = rb_simulation(model,reduced_data);
326 disp([
'time for online phase: t = ',num2str(t)]);
328 disp(
'full simulation:')
330 sim_data = detailed_simulation(model,model_data);
332 disp(['time for detailed simulation: t = ',num2str(t)]);
334 case 7 % training-error landscape
335 disp('to be adjusted to new syntax!!!');
336 disp('warning: takes a few hours!');
337 load(fullfile(rbmatlabtemp,detailedfname));
339 % runs = {
'train_linRB',
'test_linRB'};
340 % dirnames = {testdir};
341 % testdir =
'chemnitz_test_100';
343 %runs = {
'train',
'test'};
345 %dirnames = {params.RB_detailed_train_savepath, testdir};
346 dirnames = {testdir};
349 % dirnames = {params.RB_detailed_train_savepath};
351 % dirnames = {testdir};
354 %dirnames = {testdir};
356 % generate test data
if not existent
357 if ismember(
'test',runs)
358 % generate test data
if not existent
359 rand(
'state',123456);
360 M = rand_uniform(10,model.mu_ranges);
364 % load demo_nonlin_evol_detailed_data_LE_on_RB;
365 offline_data = rb_offline_prep(detailed_data,model);
368 Nmax = size(detailed_data.RB,2);
370 Ns = round((0:(Nsamples-1)) * ...
371 (Nmax-1)/(Nsamples-1)+1);
375 Mmax = size(detailed_data.BM{1},2);
376 %Mmax = 300; % size(detailed_data.BM{1},2);
379 Ms = round((0:(Msamples-1)) * ...
380 (Mmax-1)/(Msamples-1)+1);
385 for ru = 1:length(runs);
386 disp([
'starting run ',runs{ru}]);
387 % storage
for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
388 errs = zeros(Msamples,Nsamples);
389 inds = zeros(Msamples,Nsamples);
390 mu_inds = zeros(Msamples,Nsamples);
391 % assuming, that data is in the following directory
392 settings = load(fullfile(rbmatlabtemp,...
393 dirnames{ru},
'settings.mat'));
395 for mu_ind = 1:size(settings.M,2);
396 disp([
'processing parameter vector ',num2str(mu_ind),...
397 '/',num2str(size(settings.M,2))]);
398 model = set_mu(settings.M(:,mu_ind),model);
401 for Nind = 1:Nsamples
403 for Mind = 1:Msamples
407 % nonlinear simulation
410 reduced_data = rb_online_prep(offline_data,model);
412 simulation_data = rb_simulation(reduced_data,model);
413 U = rb_reconstruction(detailed_data,simulation_data);
415 [err,ind] = max(l2_errors);
416 if err(1)>errs(Mind,Nind)
417 errs(Mind,Nind) = err(1);
418 inds(Mind,Nind) = ind(1);
419 mu_inds(Mind,Nind) = mu_ind;
422 errs(Mind,Nind) = NaN;
423 inds(Mind,Nind) = NaN;
424 mu_inds(Mind,Nind) = mu_ind;
429 disp('saving temporary workspace');
430 save(fullfile(rbmatlabtemp,...
431 [results_fname_base,'_',runs{ru},
'_tmp.mat']));
433 % define stability-region as error being smaller than
434 % sqrt(diffmax * area), e.g. diffmax = 4, area = 2e-7
437 C = ones(size(errs));
438 i = find(errs<stab_limit);
445 if isempty(find(C==1)); % i.e. all stable
447 else % some not stable
448 surf(Ms, Ns,errs',C
');
459 % figure, pcolor(Ms, Ns,C);
460 title([runs{ru},' L-infty([0,T],L2) error
']);
461 set ( gca, [error_axis, 'scale
'], 'log
' );
464 eval([n_axis,'label(
''N
'');
']);
465 eval([error_axis,'label(
''error
'');
']);
473 title([runs{ru},
' time index of maximum l2-error']);
474 %set(gca,
'Zscale',
'log');
477 eval([n_axis,
'label(''N'');']);
478 %zlabel(
'time index');
479 eval([error_axis,
'label(''time index'');']);
481 save(fullfile(rbmatlabtemp,...
482 [results_fname_base,
'_',runs{ru},
'.mat']));
485 case 8 % generate figures and runtime table
487 % detailed_simulation: shock to right, v = 0.9
491 model.c_dir_left = 0;
492 model.c_dir_right = 0.5;
494 U1 = detailed_simulation(detailed_data.grid,model);
496 model.c_dir_left = 0.5;
497 model.c_dir_right = 1;
498 model.flux_vx = -0.5;
499 U2 = detailed_simulation(detailed_data.grid,model);
503 plot_element_data(detailed_data.grid,U1(:,1),model);
504 title(
'Initial Data');
505 saveas(gcf,
'shock_initial_data1',
'epsc')
507 plot_element_data(detailed_data.grid,U2(:,1),model);
508 title('Initial Data');
509 saveas(gcf,'shock_initial_data2','epsc')
511 plot_element_data(detailed_data.grid,U1(:,end),model);
513 saveas(gcf,'shock_end_data1','epsc')
515 plot_element_data(detailed_data.grid,U2(:,end),model);
517 saveas(gcf,'shock_end_data2','epsc')
519 % interpolation points
522 model.M = length(detailed_data.TM);
523 rb_plot_interpolation_points(detailed_data,model);
525 cmap = cmap(1:150,:); % take dark shades
526 cmap(1,:) = [1.0,1.0,1.0]; % light yellow
531 saveas(gcf,'shock_interpolation_points','epsc')
533 offline_data = rb_offline_prep(detailed_data,model);
534 figure, plot(offline_data.grid_local_ext);
538 set(gca,'YLim',[0,1]);
539 saveas(gcf,'shock_subgrid','epsc')
541 % figure of error convergence
543 % load(fullfile(rbmatlabtemp,detailedfname));
544 load(fullfile(rbmatlabtemp,...
545 [results_fname_base,'_test.mat']));
552 % figure, pcolor(Ms, Ns,C');
553 title([runs{ru},
' L-infty([0,T],L2) error']);
554 set( gca, [error_axis,
'scale'],
'log' );
557 eval([n_axis,
'label(''N'');']);
558 eval([error_axis,
'label(''error'');']);
559 set(gcf,
'Position',[ 120 524 496 370]);
561 tmp = load('redgreencolmap.mat');
563 % cmap = cmap(end:-1:1,:);
565 % disp('adjust colorbar!');
566 set(gcf,'Color',[1.0,1.0,1.0]);
567 saveas(gcf,'shock_error_convergence','epsc')
569 % f1 = figure, surf(Ms, Ns,errs');
570 % title(['test L^\infty([0,T],L^2) error'],'Fontsize',15);
571 % title(['test error'],'Fontsize',15); %L^\infty([0,T],L^2) error'],'Fontsize',15);
572 % set(gca,'Zscale','log');
573 % xlabel('M','Fontsize',15);
574 % ylabel('N','Fontsize',15);
575 % zlabel('error','Fontsize',15);
579 %offline_data = rb_offline_prep(detailed_data,model);
580 %plot(offline_data.grid_local_ext);
584 %set(gcf,'Color',[1,1,1]);
588 % generate runtime table
590 % load(fullfile(rbmatlabtemp,detailedfname));
591 load(fullfile(rbmatlabtemp,...
592 [results_fname_base,'_test.mat']));
596 Ns = [10,20,30,40,50,60, 70, 80, 90, 100];
597 % Ms = [15,30,45,60,75,90,105,120, 135,150];
598 Ms = 100 * ones(size(Ns));
599 t_detailed = zeros(nreps,1);
601 % model.flux_pdeg = rand(1)+1;
602 mu = rand_uniform(1,model.mu_ranges);
603 model = set_mu(mu,model);
604 % offline_data = rb_offline_prep(detailed_data,model);
606 u = detailed_simulation(model);
607 t_detailed(rep) = toc;
608 disp(['runtime for detailed: ',num2str(t_detailed(rep)),' sec.']);
611 t_reduced = zeros(nreps,length(Ns));
612 offline_data = rb_offline_prep(detailed_data,model);
613 for rn = 1:length(Ns)
616 reduced_data = rb_online_prep(offline_data,model);
618 mu = rand_uniform(1,model.mu_ranges);
619 model = set_mu(mu,model);
621 simulation_data = rb_simulation(reduced_data,model);
622 t_reduced(rep,rn) = toc;
623 disp(['N= ',num2str(model.N),', M= ',num2str(model.M)]);
624 disp(['runtime for reduced: ',num2str(t_reduced(rep,rn)),' sec.']);
628 disp(' t_mean t_std ');
629 disp([num2str(mean(t_detailed)),' ',num2str(std(t_detailed))]);
632 disp('N M t_mean t_std');
633 o = [Ns(:), Ms(:), mean(t_reduced,1)', std(t_reduced,1)'];
634 s = sprintf('%2.0f %2.0f %2.2f %2.2f \n ',o');
638 % nicely arranged error landscape figures:
640 % load(fullfile(rbmatlabtemp,...
641 % [results_fname_base,'_train.mat']));
642 load(fullfile(rbmatlabtemp,...
643 [results_fname_base,'_test.mat']));
648 load(fullfile(rbmatlabtemp,detailedfname));
649 % model = riemann_burgers_model;
650 % model.show_colorbar = 1;
651 model.c_dir_left = -1;
652 model.c_dir_right = 0;
657 model.enable_error_estimator = 1;
660 case 10 % start
demo_rb_gui without error estimation
661 load(fullfile(rbmatlabtemp,detailedfname));
662 % model = riemann_burgers_model;
663 % model.show_colorbar = 1;
667 model.enable_error_estimator = 0;
671 error('step-number is unknown!');
function [ sim_data , tictoc ] = load_detailed_simulation(m, savepath, params)
load single trajectory of previously saved results.
function Udirichlet = dirichlet_values(model, X, Y)
UDIRICHLET = DIRICHLET_VALUES([X],[Y], MODEL) Examples dirichlet_values([0,1,2],[1,1,1],struct(name_dirichlet_values, homogeneous, ... c_dir, 1)) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, xstripes, ... c_dir, [0 1 2], ... dir_borders, [0.3 0.6])) dirichlet_values([0:0.1:1],[0],struct(name_dirichlet_values, box, ... c_dir, 1, ... dir_box_xrange, [0.3 0.6], ... dir_box_yrange, [-0.1 0.1]))
function demo_rb_gui(varargin)
reduced basis demo with sliders
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 l2_error = fv_l2_error(U1, U2, W)
function computing the l2-error between the two fv-functions or function sequences in U1...
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 detailed_data = rb_basis_generation(model, detailed_data)
reduced basis construction with different methods