rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
richards_fv_model.m
1 function model = richards_fv_model(params)
2 % function model = richards_fv_model(params)
3 %
4 % optional fields of params:
5 % model_size: - 'big' simulations are done on a larger grid with very
6 % small time steps
7 % - 'small' simulations are done on a smaller grid with
8 % bigger time steps
9 % model_type:
10 % - 'linear_heat_trapezoidal' on transformed domain linear heat equation
11 % with a trapezoidal geometry parametrisation. on reference
12 % domain a linear convection-diffusion-reaction problem. It
13 % still depends on empirical interpolation, however, because
14 % of non-affine parameter dependencies of the operators
15 % introduced by geometry parametrisation. Note, that
16 % diffusion term is discretized explicitly.
17 % - 'linear_heat_polynomial' on transformed domain linear heat equation
18 % with polynomial geometry parametrization. On reference
19 % domain a linear convection-diffusion-reaction problem.
20 % Note, that diffusion term is discretized explicitly.
21 % - 'nonlinear_richards_trapezoidal' on transformed domain real Richards
22 % equation example with trapezoidal geometry parametrisation
23 % and a gravity induced convection term. On reference
24 % domain a non-linear convection-diffusion-reaction problem.
25 % Note, that diffusion term is discretized explicitly.
26 % - 'linear_heat_polynomial_implicit' On transformed domain
27 % lienar heat equation with polynomial geometry
28 % parametrisation. On reference domain a linear
29 % convection-diffusion-reaction problem with non-affine
30 % parameter dependence (resolved by EI). The diffusion term
31 % is discretized implicitly.
32 % - 'test' very small heat equation example suitable for
33 % quick test runs
34 
35  model = nonlin_evol_model_default;
36  model.name = 'richards_fv_model';
37  model.rb_problem_type = 'nonlin_evol';
38  %% data that is definitely used outside of the detailed simulations
39  % mu_names = {'c_init','hill_height'};
40  % mu_ranges = {[-1,0],[-0.1,0.1]};
41  model.mu_names = {'c_dir_low','hill_height'};
42  model.mu_ranges = {[0,0.5],[0.0,0.4]};
43 
44 
45  %% data that might be of interest outside of the detailed simulations
46  model.T = 0.9;
47  model.nt = 300;
48 
49  %% finite volume settings
50  model.diffusivity_ptr = @diffusivity_homogeneous;
51  model.diffusivity_tensor_ptr = @diffusivity_tensor_richards;
52  % diffusivity_ptr = @diffusivity_linear_gradient;
53  % diff_left = 0.001;
54  % diff_right = 0.002;
55 
56 
57  model.num_conv_flux_ptr = @fv_num_conv_flux_engquist_osher;
58  model.num_diff_flux_ptr = @fv_num_diff_flux_gradient_tensor;
59  model.flux_linear = true;
60  model.name_flux = 'richards';
61 
62 
63  model.operators_diff_implicit = @fv_operators_diff_implicit_gradient_tensor;
64  model.operators_conv_implicit = @fv_operators_zero;
65  model.operators_neumann_implicit = @fv_operators_zero;
66  model.fv_expl_conv_weight = 1.0;
67  model.fv_expl_diff_weight = 1.0;
68  model.fv_expl_react_weight = 1.0;
69 
70  % lxf_lambda = 1.25;
71  % lxf_lambda = fv_search_max_lxf_lambda([],
72 
73 
74  %% define the grid
75  model = unitcube(model);
76  model.xnumintervals = 100;
77  model.ynumintervals = 80;
78 
79  % pointer to function in rbmatlab/datafunc/dirichlet_values
80  model.dirichlet_values_ptr = @dirichlet_values_uplow;
81  model.c_dir = 0.1;
82  model.dir_box_xrange = [-1 2];
83  model.dir_box_yrange = 0.5 + [1/320 3/320];
84  model.c_dir_left = 1;
85  model.c_dir_right = 1;
86  model.c_dir_up = 0;
87  model.c_dir_low = 0.5;
88  model.dir_middle = 0.7;
89  % name_dirichlet_values = 'homogeneous';
90  % c_dir = 0;
91 
92  % pointer to function in rbmatlab/datafunc/neumann_values
93  %neumann_values_ptr = @neumann_values_homogeneous;
94  model.neumann_values_ptr = @neumann_values_zero;
95  model.c_neu = 0;
96 
97  % name of function in rbmatlab/datafunc/init_values/
98  model.init_values_ptr = @init_values_transformed_blobs;
99  % parameters for data functions
100  model.blob_height = 0.1;
101 
102  model.filecache_velocity_matrixfile_extract = 0;
103 
104  model.geometry_spline_type = 'cubic';
105  model.hill_height = 0.0;
106  model.geometry_transformation_spline_x = [ 0 0.5 1 ];
107  model.geometry_transformation_spline_y = [ 0 -0.033 0 ];
108 
109  model.k = 0.005;
110  model.gravity = 0;
111 
112 
113  % settings for CRB generation
114  model.CRB_generation_mode = 'param-time-space-grid';
115  model.Mmax = 150;
116  model.sMmax = {50, 50};
117  model.M = 40;
118  model.Mstrich = 5;
119  model.ei_stop_on_Mmax = 1;
120  % build ei-basis with WHOLE trajectory
121  model.ei_target_error = 'interpol';
122  model.ei_numintervals = [2,4];
123 
124  model.RB_stop_Nmax = 50;
125  model.RB_stop_epsilon = 1e-5;
126  model.RB_stop_max_val_train_ratio = inf;
127  model.RB_numintervals = [2,4];
128 
129  % target accuracy epsilon:
130 
131  % decide, whether estimator or true error is error-indicator for greedy
132  model.RB_error_indicator = 'error'; % true error
133  % RB_error_indicator = 'estimator'; % Delta from rb_simulation
134  % RB_error_indicator = 'ei_estimator_test'; % Delta from rb_simulation testet against true error
135 
136  model.divclean_mode = false;
137  model.flux_quad_degree = 1;
138 
139  model.conv_flux_ptr = @conv_flux_linear;
140 
141  model.velocity_ptr = @velocity_richards;
142 
143  model.conv_flux_derivative_ptr = @(glob, U, params) params.velocity_ptr(glob, params);
144 
145 
146  model.data_const_in_time = 0;
147 
148  % set all to dirichlet-boundary by specifying "rectangles", all
149  % boundary within is set to boundary type by bnd_rect_index
150  model.bnd_rect_corner1 = [-0.999999, 0; ...
151  0, -0.9999999 ];
152  model.bnd_rect_corner2 = [ 2, 0.999999; ...
153  0.999999, 2 ];
154  % -1 means dirichlet, -2 means neumann
155  model.bnd_rect_index = [ -1, -2 ];
156 
157  if isfield(params, 'separate_CRBs') && params.separate_CRBs
158  model.separate_CRBs = params.separate_CRBs;
159  savepath_infix = '';
160  else
161  savepath_infix = '_one_CRB';
162  end
163 
164  model.ei_detailed_savepath = [model.name,'_',params.model_type,...
165  savepath_infix,'_ei_data_interpol'];
166  model.ei_operator_savepath = [model.name,'_',params.model_type,...
167  savepath_infix,'_ei_operators_interpol'];
168 
169  model.ei_space_operators = { model.L_E_local_ptr };
170 
171  model.geometry_transformation = 'spline';
172 
173  model.stencil_mode = 'vertex';
174  model.local_stencil_size = 2;
175 
176  model.RB_detailed_train_savepath = model.ei_detailed_savepath;
177 
178  if isfield(params, 'model_size')
179  if strcmp(params.model_size, 'big') == 1
180  model.xnumintervals = 200;
181  model.ynumintervals = 200;
182  end
183  if strcmp(params.model_size, 'small') == 1
184  model.xnumintervals = 50;
185  model.ynumintervals = 50;
186  model.T = 2;
187  model.nt = 20;
188  end
189  end
190  if isfield(params, 'model_type')
191  if strcmp(params.model_type, 'nonlinear') == 1
192  model.geometry_spline_type = 'cubic';
193  model.nt = 300;
194  elseif strcmp(params.model_type, 'linear_heat_trapezoidal') == 1
195  model.geometry_spline_type = 'affine';
196  model.T = 0.5;
197  model.nt = 260;
198  model.c_init = 0;
199  elseif strcmp(params.model_type, 'richards_affine') == 1
200  model.geometry_spline_type = 'affine';
201  model.diffusivity_ptr = @diffusivity_richards_nonlinear;
202  model.init_values_ptr = @init_values_transformed_blobs_richards;
203  % set parameters for nonlinear gradient
204  model.k = 0.002;
205  model.gravity = 0.5;
206  model.clim = [0.22,0.40];
207  % set new time parameters
208  model.T = 0.01;
209  model.nt = 195;
210  model.c_dir_up = 0.32;
211  model.c_dir_low = 0.35;
212  model.xnumintervals = 80;
213  model.ynumintervals = 100;
214  % set default mu values
215  model.c_init = 0.34;
216  model.hill_height = 0.50;
217  % choose different mu vector
218  model.mu_names = {'hill_height','c_init'};
219  model.mu_ranges = {[0,0.5],[0.22,0.34]};
220  % only neumann boundaries
221  model.bnd_rect_index = [ -2, -2 ];
222 
223  model.ei_numintervals = [5,6];
224  model.RB_numintervals = [5,6];
225  ht = @(t)(t-0.218)./(0.52-0.218);
226  model.richards_perm_ptr = @(t)(2.1*ht(t).^0.5.*(1-(1-ht(t).^(7/6)).^(6/7)).^2);
227  model.richards_retention_ptr = @(t)(35.2574./(1./(3.311258.*t-0.7218543).^(6/7)-1).^(6/7)./(3.311258.*t-0.721854).^(13/7));
228  model.postprocess = @postprocess_gravity;
229  model.RB_stop_Nmax = 10;
230  elseif strcmp(params.model_type, 'implicit_nonaffine_linear') == 1
231  model.implicit_nonlinear = true;
232  model.hill_height = 0.4;
233  model.geometry_spline_type = 'cubic';
234  % model.mu_names = {'hill_height', 'blob_height'};
235  model.mu_ranges = {[0,0.4], [0.0,0.3]};
236  % model.c_init = 0;
237  % model.bnd_rect_index = [ -2, -2 ];
238  model.k = 0.005;
239  % k = 0.005;
240  model.nt = 100;
241  model.fv_expl_diff_weight = 0.0;
242  model.fv_impl_diff_weight = 1.0;
243  model.data_const_in_time = 1;
244  model.T = 1.5;
245  model.ei_space_operators = { model.L_E_local_ptr, model.L_I_local_ptr };
246  elseif strcmp(params.model_type, 'test')
247  model.T = 0.15;
248  model.nt = 30;
249  model.Mmax = 20;
250  model.RB_stop_Nmax = 5;
251  model.M = model.Mmax;
252  model.Mstrich = 0;
253  model.geometry_spline_type = 'affine';
254  model.data_const_in_time = true;
255  model.N = model.RB_stop_Nmax;
256  model.ei_numintervals = [2,2];
257  model.RB_numintervals = [2,2];
258  else
259  error(['selected model type "', params.model_type, '" is unknown.']);
260  end
261  end
262  model.model_type = params.model_type;
263 
264  model = model_default(model);
265 
266  model.verbose = 1;
267  model.debug = 0;
268 
269 end
270 % vim: set et sw=2:
271 %| \docupdate
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 data = postprocess_gravity(data, glob, params)
subtracts a previously added addent induced by gravitational effects.
function num_flux = fv_num_diff_flux_gradient_tensor(model, model_data, U, NU_ind)
computes a numerical diffusive flux for a diffusion problem including a tensor
function model = model_default(model, T, nt)
model = model_default(model)
Definition: model_default.m:17
function model = unitcube(model)
function adding fields to model for generating a 2D rectgrid with 100 x 100 elements on the unit-squa...
Definition: unitcube.m:17
function num_flux = fv_num_conv_flux_engquist_osher(model, model_data, U, NU_ind)
Function computing a numerical convective Engquist-Osher flux matrix.
function [ flux , lambda ] = conv_flux_linear(glob, U, params)
function computing the convective flux of a convection problem.