1 % small script demonstrating the richards equation with
explicit fv
2 % discretization and RB model reduction
3 % Bernard Haasdonk 14.8.2007
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%%%% Select here, what is to be performed with the richards data %%%%%%%%
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %step = 1 % single detailed simulation with given data and plot. Run
9 %
this with varying parameters mu until sure that scheme
10 % is stable. Modify dt or the data-functions accordingly,
11 % until a nice parameter-domain with uniformly stable
12 % detailed scheme is obtained.
13 %step = 2 % generate colateral reduced basis of L_E operator
14 %step = 3 % compare single detailed simulation with and without
15 % interpolated space operator
16 %step = 4 % generate dummy reduced basis from single trajectory and check, if
17 % ei_interpolation with projection on this space maintains
18 % result. A simple reduced simulation can also be
19 % performed. All results should be visually identical
20 %step = 5 % generate reduced basis
21 %step = 6 % time measurement of reduced simulation and
22 % use reduced basis in rb_demo_gui
23 %step = 7 % generate error-landscape over varying N and M
24 % can take several hours!!!
25 %step = 8 % do runtime comparisons between detailed and reduced simulations
26 %step = 10 % write out some nice pictures
31 imsavepath =
'/data/dune_work/private/m_droh01/images';
33 for si=1:length(steps)
37 % params.model_type =
'implicit_nonaffine_linear';
38 params.model_type =
'linear_heat_trapezoidal';
39 %params.model_size =
'big';
40 params.model_size =
'small';
41 % params.model_size =
'small';
42 % params.model_type =
'nonlinear';
43 params.separate_CRBs =
false;
46 params.RB_detailed_test_savepath =
'test_data_100';
47 params.RB_test_size = 100;
51 params.step7_outputfile = fullfile(rbmatlabresult, [params.model_type,
'_step7_output']);
52 params.step8_outputfile = fullfile(rbmatlabresult, [params.model_type,
'_step8_output']);
53 extrapar.Mstrich = 20;
57 % email to which status messages about processing mu vectors is sent
58 % params.email =
'mdrohmann@uni-muenster.de';
60 %% parameters
for visualization
61 plot_params.show_colorbar = 1;
62 plot_params.colorbar_mode =
'EastOutside';
63 plot_params.no_lines =
true;
64 % plot_params.clim = [0, 0.5];
65 % plot_params.clim_tight =
true;
67 plot_params.bind_to_model =
true;
68 plot_params.yscale_uicontrols = 0.6;
70 % output-filenames in rbmatlabresult
71 CRBfname = ['richards_fv_', params.model_type, '_CRB_interpol.mat'];
72 detailedfname = ['richards_fv_', params.model_type, '_detailed_interpol.mat'];
75 case 0 % initialize model data
76 model = richards_fv_model(params);
77 params.mu_ranges = model.mu_ranges;
78 params.mu_names = model.mu_names;
79 mu_default = cellfun(@mean, model.mu_ranges);
81 model_data = gen_model_data(model);
82 case 1 % single detailed simulation and plot
84 case 2 % empirical interpolation of space
operator
86 case 3 % compare detailed simulation with and without interpolation
87 load(fullfile(rbmatlabresult,CRBfname));
89 case 4 % construct dummy reduced basis by single trajectory and simulate
90 tmp=load(fullfile(rbmatlabresult,CRBfname));
91 detailed_data=tmp.detailed_data;
93 case 5 % reduced basis
94 load(fullfile(rbmatlabresult,CRBfname));
97 load(fullfile(rbmatlabresult,detailedfname));
99 case 7 % training-error landscape & time comparisons
100 load(fullfile(rbmatlabresult,detailedfname));
103 % maximum number of reduced basis functions
104 model.Nmax = size(detailed_data.RB,2);
105 % maximum number of collateral reduced basis functions
106 model.Mmax = size(detailed_data.BM{1},2);
107 % number of reduced basis sizes to be tested.
109 % number of collateral reduced basis sizes to be tested.
111 % sample of values
for the coupling variable
'c' for which errors are
114 % axis description
for plots with coupling variable
'c'
115 cdescr =
'(N,M) = c * (Nmax, Mmax)';
116 % size of mu vector test
set.
122 load(params.step7_outputfile);
123 for i = 1:length(output)
127 disp(
'warning: takes a few hours!');
128 load(fullfile(rbmatlabresult,detailedfname));
130 % maximum number of reduced basis vectors
131 model.Nmax = model.RB_stop_Nmax;
132 % maximum number of collateral reduced basis vectors used
for
133 model.Mmax = size( detailed_data.BM{1}, 2 ) - extrapar.Mstrich;
135 % number of reduced basis sizes to be tested.
137 % number of collateral reduced basis sizes to be tested.
139 % Mstrich values to be tested
140 Mstrich_samples = 1:extrapar.Mstrich;
141 % sample of values
for the coupling variable
'c' for which
142 % estimators are computed.
143 csamples = 0.1:0.1:1;
144 % size of mu vector test
set
146 % axis description
for plots with coupling variable
'c'
147 cdescr =
'(N,M) = c * (Nmax, Mmax)';
151 load(params.step8_outputfile);
152 for i = 1:length(output)
153 output{i}.bound = 1500;
154 output{i}.stab_limit = 10000;
158 load(fullfile(rbmatlabresult,detailedfname));
159 model_data = detailed_data;
160 testdir =
'richards_test_10';
162 model.Nmax = size(detailed_data.RB,2);
164 % generate test data
if not existent
165 rand(
'state',123456);
166 M = rand_uniform(10,params.mu_ranges);
169 % load demo_nonlin_evol_detailed_data_LE_on_RB;
170 reduced_data = gen_reduced_data(model,detailed_data);
172 settings = load(fullfile(rbmatlabtemp,...
173 testdir,
'settings.mat'));
177 % storage
for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
178 errs = zeros(Msamples,Nsamples);
179 model.RB_error_indicator=
'ei_estimator_test';
181 Mstrich_set = 1:10:100;%[1,2,3,4,5,10,15,20];
183 deltas = zeros(size(M,2),length(Mstrich_set));
184 for mu_ind = 1:size(M,2);
185 disp([
'processing parameter vector ',num2str(mu_ind),...
186 '/',num2str(size(M,2))]);
187 for msi = 1:length(Mstrich_set)
189 model = model.set_mu(model,settings.M(:,mu_ind));
195 mtemp.M = extrapar.M;
196 mtemp.Mstrich = Mstrich_set(msi);
197 mtemp.N = extrapar.N;
198 simulation_data = rb_simulation(mtemp,reduced_data);
201 deltas(mu_ind,msi) = simulation_data.Delta(end);
202 simulation_data = rb_reconstruction(mtemp,detailed_data,simulation_data);
203 error = model.error_algorithm(simulation_data.U, sim_data.U, ...
204 detailed_data.W, mtemp);
205 disp(['true error: ', num2str(max(error))]);
207 %
plot_sim_data(model,detailed_data,simulation_data,plot_params);
210 plot(Mstrich_set,deltas(mu_ind,:));
214 load('datafiles/strange.mat');
216 tmp = load(fullfile(rbmatlabtemp,detailedfname));
217 detailed_data = tmp.detailed_data;
218 plot_params.no_lines = 1;
220 plot_params.no_axes = 1;
221 plot_params.transparent_background = 1;
222 plot_params.show_colorbar = 0;
223 %plot_params.shrink_factor = 1.13;
224 plot_params.colormap = Cmap;
225 plot_params.clim = [0.0,0.35];
227 reduced_data = gen_reduced_data(model, detailed_data);
228 params.N = 14;%size(detailed_data.RB,2);
229 params.M = 100;%size(detailed_data.QM,2);
230 reduced_data = extract_reduced_data_subset(model, reduced_data, params);
232 model.hill_height = 0;
234 cparams.range = cell(length(model.mu_ranges)+1,1);
235 cparams.range{1} = [1, model.nt];
236 cparams.range(2:end) = model.mu_ranges(:);
237 cparams.numintervals = 2*ones(length(model.mu_ranges)+1,1);
240 parammat = cgrid.vertex;
241 parammat(:,1) = round(parammat(:,1));
243 for i=1:length(parammat)
244 tstep = parammat(i,1);
245 newmodel = newmodel.set_mu(model, parammat(i,2:end));
246 rb_sim_data = rb_simulation(newmodel, reduced_data, params);
247 rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
249 % transformed version
250 plot_element_data(newmodel, model_data.grid, rb_sim_data.U(:,tstep), plot_params);
252 mustring = strrep(strrep(strrep(strrep(mat2str(parammat(i,2:end)),' ', '_'),'[',''),']',''),'.','p');
253 fp = fullfile(imsavepath, params.model_type);
257 fp = fullfile(fp, ['t_', num2str(tstep)]);
261 fn = ['sample_mu_', mustring];
262 disp(['Processing ', fn, ' ...']);
264 set(gcf,'Color','none');
265 set(gca,'Color','none');
266 set(gcf,'InvertHardCopy','off');
269 ranges{1} =
get(gca,
'XLim');
270 ranges{2} =
get(gca,
'YLim');
272 export_fig(fullfile(fp, fn),
'-pdf',
'-eps',
'-png');
275 colorbarfn = fullfile(imsavepath, params.model_type,
'colorbar');
276 if ~exist([colorbarfn,
'.pdf'],
'file')
278 yticks = get(colorh, 'YTick');
279 set(colorh, 'YTick',[]);
280 set(colorh, 'YTickLabel',[]);
282 export_fig(colorbarfn, '-pdf', '-eps', '-png', colorh);
283 cbim = export_fig(colorh);
287 tikzfile(newmodel, fp, fn, ranges, size(im), 5);
290 % untransformed version
291 fpu = fullfile(imsavepath, params.model_type, 'untransformed');
295 fpu = fullfile(fpu, ['t_', num2str(tstep)]);
300 plot_element_data(model, model_data.grid, rb_sim_data.U(:,tstep), plot_params);
302 set(gcf,'Color','none');
303 set(gca,'Color','none');
304 set(gcf,'InvertHardCopy','off');
306 ranges{1} =
get(gca,
'XLim');
307 ranges{2} =
get(gca,
'YLim');
309 export_fig(fullfile(fpu, fn),
'-pdf',
'-eps',
'-png');
311 tikzfile(newmodel, fpu, fn, ranges, size(im), 5);
316 load(
'richards_affine_detailed_interpol');
318 rdata = gen_reduced_data(model, detailed_data);
320 model.N = size(detailed_data,2);
323 rdata2 = extract_reduced_data_subset(model, rdata);
325 rdata1 = extract_reduced_data_subset(model, rdata);
327 TM1 = rdata1.TM_local{1}(1:10);
328 TM2 = rdata2.TM_local{1}(1:10);
330 gl1 = rdata1.grid_local_ext{1};
331 gl2 = rdata2.grid_local_ext{1};
333 if(any(gl1.CX(TM1) ~= gl2.CX(TM2)))
336 if(any(gl1.CY(TM1) ~= gl2.CY(TM2)))
339 if(any(any(gl1.X(gl1.VI(TM1,:)) ~= gl2.X(gl2.VI(TM2,:)))))
342 if(any(any(gl1.Y(gl1.VI(TM1,:)) ~= gl2.Y(gl2.VI(TM2,:)))))
345 if(any(any(gl1.SX(TM1,:) ~= gl2.SX(TM2,:))))
348 if(any(any(gl1.SY(TM1,:) ~= gl2.SY(TM2,:))))
351 if(any(gl1.EL(TM1,:) ~= gl2.EL(TM2,:)))
354 if(any(gl1.DC(TM1,:) ~= gl2.DC(TM2,:)))
357 if(any(gl1.NX(TM1,:) ~= gl2.NX(TM2,:)))
360 if(any(gl1.NY(TM1,:) ~= gl2.NY(TM2,:)))
363 if(any(gl1.ECX(TM1,:) ~= gl2.ECX(TM2,:)))
366 if(any(gl1.ECY(TM1,:) ~= gl2.ECY(TM2,:)))
370 NBI1 = gl1.NBI(TM1,:);
371 NBI2 = gl2.NBI(TM2,:);
372 if(any(gl1.CX(NBI1) ~= gl2.CX(NBI2)))
375 if(any(gl1.CY(NBI1) ~= gl2.CY(NBI2)))
378 if(any(any(gl1.X(gl1.VI(NBI1,:)) ~= gl2.X(gl2.VI(NBI2,:)))))
381 if(any(any(gl1.Y(gl1.VI(NBI1,:)) ~= gl2.Y(gl2.VI(NBI2,:)))))
384 if(any(any(gl1.SX(NBI1,:) ~= gl2.SX(NBI2,:))))
387 if(any(any(gl1.SY(NBI1,:) ~= gl2.SY(NBI2,:))))
390 for i = 1:length(TM1)
394 error('step-number is unknown!');