1 % small script demonstrating the burgers equation with
explicit fv
2 % discretization and RB model reduction
5 % Bernard Haasdonk 14.8.2007
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %%%%%% Select here, what is to be performed with the burgers data %%%%%%%%%%
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10 %step = 1 % single detailed simulation with given data and plot. Run
11 %
this with varying parameters mu until sure that scheme
12 % is stable. Modify dt or the data-functions accordingly,
13 % until a nice parameter-domain with uniformly stable
14 % detailed scheme is obtained.
15 %step = 2 % generate colateral reduced basis of L_E operator
16 step = 3 % compare single detailed simulation with and without
17 % interpolated space operator
18 %step = 4 % generate dummy reduced basis from single trajectory and check, if
19 % ei_interpolation with projection on this space maintains
20 % result. A simple reduced simulation can also be
21 % performed. All results should be visually identical
22 %step = 5 % generate reduced basis
23 %step = 6 % time measurement of reduced simulation and
24 % use reduced basis in rb_demo_gui
25 %step = 7 % generate error-landscape over varying N and M
26 % can take several hours!!!
28 % output-filenames in rbmatlabresult
29 CRBfname = 'burgers_fv_CRB_interpol.mat';
30 detailedfname =
'burgers_fv_detailed_interpol.mat';
34 % rectangular or triangular
35 params.gridtype =
'triagrid';
36 params.grid_initfile =
'2dsquaretriang';
38 %params.gridtype =
'rectgrid';
39 %params.xnumintervals = 20;
40 %params.ynumintervals = 20;
41 %params.xrange = [0,1];
42 %params.yrange = [0,1];
44 %
set all to neuman-boundary by specifying "rectangles", all
45 % boundary within is set to boundary type
46 params.bnd_rect_corner1 = [-1; ...
48 params.bnd_rect_corner2 = [ 2; ...
51 % -1 means dirichlet, -2 means neumann
52 params.bnd_rect_index = [-2]; % first is dirichlet, second neuman
56 %params.nt = 200; % guess and adjust later
59 params.nt = 80; % guess and adjust later
60 %params.nt = 100; % guess and adjust later
62 %params.nt = 50; % guess and adjust later
63 %params.400 => non-bounded u
67 params.axis_equal = 1;
71 %params.name_init_values =
'blobs';
74 %params.beta = 1; % init data weight between 0 and 1
76 % name of
function in rbmatlab/datafunc/init_values
77 params.name_init_values =
'gradient_box';
78 % parameters
for data functions
79 params.c_init_up = 1.0;
80 %params.c_init_lo = 1.0;
81 params.c_init_lo = 0.0;
82 %params.c_init_lo = -1.0;
85 %params.name_dirichlet_values =
'homogeneous';
88 % name of
function in RBmatlab/datafunc/neuman_values
89 params.name_neuman_values =
'homogeneous';
92 %params.name_neuman_values =
'homogeneous';
96 params.name_flux =
'burgers_parabola';
97 %
for simplifying computation in
case of linear flux:
98 params.flux_linear = 0;
99 %params.flux_quad_degree = 2;
101 % parameter of chosen conv-flux
103 %params.vrot_angle = 0;
104 params.vrot_angle = -pi/4;
105 %params.vrot_angle = -pi/10;
106 params.flux_pdeg = 2; % burgers exponent
109 params.mu_names = {
'c_init_lo',
'vrot_angle'};
110 params.mu_ranges = {[-1,1],[-pi/4, 0]};
111 params.rb_problem_type =
'nonlin_evol';
112 %params.problem_type =
'lin_evol';
114 % finite volume settings
115 params.name_diffusive_num_flux =
'none';
116 params.name_convective_num_flux =
'lax-friedrichs';
117 params.lxf_lambda = 1.66
118 %params.lxf_lambda = 1.25;
119 %params.lxf_lambda = fv_search_max_lxf_lambda([],params)
120 %params.name_convective_num_flux =
'enquist-osher';
122 % projection of analytical initial data to discrete
function
123 params.init_values_algorithm =
'fv_init_values';
124 params.implicit_operators_algorithm =
'fv_operators_implicit';
125 %
get mass matrix
for inner product computation from DOF-vectors
126 params.inner_product_matrix_algorithm = 'fv_inner_product_matrix';
127 params.data_const_in_time = 1;
129 % settings
for empirical interpolation
130 params.L_E_local_name =
'fv_conv_explicit_space';
131 par.range = params.mu_ranges;
132 params.ei_numintervals = [4,4]; % 8x8 parameter grid
135 params.ei_detailed_savepath =
'burgers_fv_5x5_detailed';
136 params.ei_operator_savepath =
'burgers_fv_5x5_ei_operator';
137 params.ei_time_indices = 1:params.nt;
138 params.CRB_generation_mode =
'param-time-space-grid';
139 params.ei_target_error =
'interpol';
141 % settings
for reduced basis generation
142 params.detailed_simulation_algorithm =
'detailed_simulation';
143 %params.operators_algorithm =
'fv_operators_implicit_explicit';
144 %params.init_values_algorithm =
'fv_init_values';
145 %params.inner_product_matrix_algorithm =
'fv_inner_product_matrix';
146 params.error_algorithm =
'fv_error';
147 %params.lxf_lambda = 1.0194e+003;
148 %params.data_const_in_time = 1;
149 params.error_norm =
'l2';
150 params.RB_extension_algorithm =
'RB_extension_PCA_fixspace';
151 %
'RB_extension_max_error_snapshot'};
152 params.RB_stop_timeout = 2*60*60; % 2 hours
153 params.RB_stop_epsilon = 1e-5;
154 params.RB_error_indicator =
'error'; %
true error
155 %params.RB_error_indicator =
'estimator'; % Delta from rb_simulation
156 params.RB_stop_Nmax = 100;
157 params.RB_generation_mode =
'uniform_fixed';
158 params.RB_numintervals = params.ei_numintervals;
159 params.RB_detailed_train_savepath = params.ei_detailed_savepath;
162 case 1 % single detailed simulation and plot
163 disp(
'performing single detailed simulation')
164 grid = construct_grid(params);
165 U = detailed_simulation(grid, params);
166 plot_element_data_sequence(grid,U,params);
167 case 2 % empirical interpolation of space operator
168 disp('constructing colateral reduced basis')
169 params.RB_generation_mode = 'none';
170 detailed_data = rb_detailed_prep(params);
171 save(fullfile(rbmatlabresult,CRBfname),...
172 'detailed_data','params');
173 detailed_data.RB = zeros(detailed_data.grid.nelements,1);
174 rb_plot_detailed_data(detailed_data,params);
176 case 3 % compare detailed simulation with and without interpolation
177 disp(['comparison between detailed simulation with and without', ...
179 tmp = load(fullfile(rbmatlabresult,CRBfname));
180 detailed_data = tmp.detailed_data;
181 U = detailed_simulation(detailed_data.grid, params);
182 params.title = 'detailed simulation without interpolation';
183 plot_element_data_sequence(detailed_data.grid,U,params);
184 params.M = size(detailed_data.QM{1},2);
186 Ui = detailed_ei_nonlin_evol_simulation(detailed_data,params);
187 params.title =
'detailed simulation with empirical interpolation';
188 plot_element_data_sequence(detailed_data.grid,Ui,params);
189 params.title =
'error';
190 plot_element_data_sequence(detailed_data.grid,abs(Ui-U),params);
191 disp([
'maximum l-infty error:',num2str(max(max(abs(Ui-U))))]);
192 case 4 % construct dummy reduced basis by single trajectory and simulate
193 load(fullfile(rbmatlabresult,CRBfname));
194 disp(
'detailed interpolated simulation for basis construction:')
195 params.M = size(detailed_data.QM{1},2);
196 Ui = detailed_ei_nonlin_evol_simulation(detailed_data,params);
197 A = feval(params.inner_product_matrix_algorithm,detailed_data.grid,params);
198 UON = orthonormalize(Ui,A);
199 detailed_data.RB = UON;
200 disp(
'reduced simulation:')
201 offline_data = rb_offline_prep(detailed_data,params);
202 params.N = size(detailed_data.RB,2);
203 params.M = size(detailed_data.QM{1},2);
204 reduced_data = rb_online_prep(offline_data, params);
205 simulation_data = rb_simulation(reduced_data,params);
206 Uappr = rb_reconstruction(detailed_data,simulation_data);
208 disp(
'detailed interpolated and rb-projected simulation:')
209 params.M = size(detailed_data.QM{1},2);
210 params.N = size(detailed_data.RB,2);
211 Uirb = detailed_ei_rb_proj_nonlin_evol_simulation(detailed_data,params);
213 params.title =
'reduced simulation result';
214 plot_element_data_sequence(detailed_data.grid,Uappr,params);
215 params.title =
'detailed ei_interpol simulation result';
216 plot_element_data_sequence(detailed_data.grid,Ui,params);
217 params.title =
'detailed ei_interpol and rb_projected simulation result';
218 plot_element_data_sequence(detailed_data.grid,Uirb,params);
220 case 5 % reduced basis
221 disp(
'constructing reduced basis')
222 tmp = load(fullfile(rbmatlabresult,CRBfname));
223 detailed_data = tmp.detailed_data;
226 save(fullfile(rbmatlabresult,detailedfname),...
227 'detailed_data','params');
228 plot(detailed_data.RB_info.max_err_sequence);
229 set(gca,'Yscale','log');
230 title('RB-generation error convergence');
232 load(fullfile(rbmatlabresult,detailedfname));
233 disp('reduced simulation:')
234 offline_data = rb_offline_prep(detailed_data,params);
235 params.N = size(detailed_data.RB,2);
236 params.M = size(detailed_data.QM{1},2);
237 reduced_data = rb_online_prep(offline_data, params);
239 simulation_data = rb_simulation(reduced_data,params);
241 disp([
'time for online phase: t = ',num2str(t)]);
243 disp(
'full simulation:')
245 U = detailed_simulation(detailed_data.grid, params);
247 disp(['time for detailed simulation: t = ',num2str(t)]);
251 case 7 % training-error landscape
252 disp('warning: takes a few hours!');
253 load(fullfile(rbmatlabresult,detailedfname));
255 % runs = {
'train_linRB',
'test_linRB'};
256 % dirnames = {testdir};
257 % testdir =
'chemnitz_test_100';
259 % runs = {
'train',
'test'};
260 % dirnames = {params.RB_detailed_train_savepath, testdir};
263 dirnames = {params.RB_detailed_train_savepath};
266 %dirnames = {testdir};
268 % generate test data
if not existent
269 if ismember(
'test',runs)
270 % generate test data
if not existent
271 rand(
'state',123456);
272 M = rand_uniform(100,params.mu_ranges);
276 % load demo_nonlin_evol_detailed_data_LE_on_RB;
277 offline_data = rb_offline_prep(detailed_data,params);
280 Nmax = size(detailed_data.RB,2);
282 Ns = round((0:(Nsamples-1)) * ...
283 (Nmax-1)/(Nsamples-1)+1);
287 Mmax = size(detailed_data.BM,2);
288 %Mmax = 300; % size(detailed_data.BM,2);
289 Ms = round((0:(Msamples-1)) * ...
290 (Mmax-1)/(Msamples-1)+1);
292 for ru = 1:length(runs);
293 disp(['starting run ',runs{ru}]);
294 % storage
for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
295 errs = zeros(Msamples,Nsamples);
296 inds = zeros(Msamples,Nsamples);
297 mu_inds = zeros(Msamples,Nsamples);
298 % assuming, that training data is in the following directory
299 settings = load(fullfile(rbmatlabtemp,...
300 dirnames{ru},
'settings.mat'));
302 for mu_ind = 1:size(settings.M,2);
303 disp([
'processing parameter vector ',num2str(mu_ind),...
304 '/',num2str(size(settings.M,2))]);
305 params = set_mu(settings.M(:,mu_ind),params);
308 for Nind = 1:Nsamples
310 for Mind = 1:Msamples
314 % nonlinear simulation
317 reduced_data = rb_online_prep(offline_data,params);
319 simulation_data = rb_simulation(reduced_data,params);
320 U = rb_reconstruction(detailed_data,simulation_data);
322 [err,ind] = max(l2_errors);
323 if err(1)>errs(Mind,Nind)
324 errs(Mind,Nind) = err(1);
325 inds(Mind,Nind) = ind(1);
326 mu_inds(Mind,Nind) = mu_ind;
329 errs(Mind,Nind) = NaN;
330 inds(Mind,Nind) = NaN;
331 mu_inds(Mind,Nind) = mu_ind;
336 disp('saving temporary workspace in burgers_fv_step7_*_tmp.mat');
337 save(fullfile(rbmatlabtemp,...
338 ['burgers_fv_step7_result_',...
339 runs{ru},
'_tmp.mat']));
341 % define stability-region as error being smaller than
342 % sqrt(diffmax * area), e.g. diffmax = 4, area = 2e-7
345 C = ones(size(errs));
346 i = find(errs<stab_limit);
349 if isempty(find(C==1)); % i.e. all stable
350 f1 = figure, surf(Ms, Ns,errs
');
351 else % some not stable
352 f1 = figure, surf(Ms, Ns,errs',C
');
355 % figure, pcolor(Ms, Ns,C);
356 title([runs{ru},' L-infty([0,T],L2) error
']);
357 set(gca,'Zscale
','log
');
362 f2 = figure, surf(Ms, Ns,inds');
363 title([runs{ru},
' time index of maximum l2-error']);
364 %
set(gca,
'Zscale',
'log');
367 zlabel(
'time index');
369 save(fullfile(rbmatlabresult,...
370 [
'burgers_fv_step7_result_',...
375 error(
'step-number is unknown!');
379 % TO BE ADJUSTED TO NEW SYNTAX