1 function [sim_data, fl] = nonlin_evol_detailed_simulation(model,model_data)
2 %
function sim_data = nonlin_evol_detailed_simulation(model,model_data)
4 % required fields of model:
6 % nt : number of time-intervals until T, i.e. nt+1
7 % solution slices are computed
8 % init_values_algorithm : name of
function for computing the
9 % initvalues-DOF with arguments (model, model_data)
10 % example: fv_init_values
11 % data_const_in_time :
if this optional field is 1, the time
12 % evolution is performed with constant operators,
13 % i.e. only the initial-time-operators are computed
14 % and used throughout the time simulation.
15 % L_E_local_ptr :
function pointer to the local space-discretization
16 %
operator evaluation with syntax
18 % INC_local = L_local(U_local_ext, ind_local,
20 % where *_local indicates results on a set of
21 % elements, *_local_ext includes information including
22 % their neighbours, i.e. the grid must be known
23 % also
for the neighbours, the values of the
24 % previous timstep must be known on the
25 % neighbours. The subset of elements, which are
26 % to be computed are given in ind_local and the
27 % values produced are INC_local.
28 % A time-step can then be realized by
29 % NU_local = U_local_ext(ind_local) - dt * INC_local
32 % Martin Drohmann 9.12.2007 based on detailed_nonlin_evol_simulation.m by
36 error(
'wrong number of parameters!');
40 disp(
'entered detailed simulation ');
43 if isempty(model_data)
44 model_data = gen_model_data(model);
47 % for a detailed simulation we do not need the affine parameter
49 model.decomp_mode = 0;
51 % initial values by midpoint evaluation
52 U0 = model.init_values_algorithm(model,model_data);
54 U = zeros(length(U0(:)),model.nt+1);
57 if ~isfield(model,'data_const_in_time')
58 model.data_const_in_time = 0;
61 if ~isfield(model, 'newton_regularisation')
62 model.newton_regularisation = 0;
65 fl = cell(model.nt, 1);
67 if model.fv_impl_diff_weight ~= 0 && ~model.newton_solver
68 % diffusive part is computed implicitly.
70 [L_diff_I, const_diff_I] = model.implicit_operators_algorithm(model, model_data, []);
74 model.dt = model.T / model.nt;
80 model.t = (t-1)*model.dt;
83 disp(['entered time-loop step ',num2str(t), ' for non linear evolution equation']);
88 [inc, fluxes] = model.L_E_local_ptr(model, model_data, U(:,t), []);
93 if all([model.fv_impl_conv_weight, model.fv_impl_diff_weight, ...
94 model.fv_impl_react_weight] == [0, 0, 0])
95 % purely
explicit scheme
96 U(:,t+1) = U(:,t) - model.dt * inc;
97 elseif model.newton_solver
98 V = zeros(length(U0(:)),model.newton_steps+1);
99 exprhs = model.dt * inc;
100 V(:,1) = U(:,t) - exprhs;
105 for n=1:model.newton_steps;
106 if model.verbose > 10
107 disp([
'Computing time step no. ', num2str(t), ...
108 ', Newton step no. ', num2str(n), ...
109 '... Residual norm: ']);
111 model.t = model.t + model.dt;
112 Urhs = model.L_I_local_ptr(model, model_data, V(:,n), []);
113 [gradL_I, gL_I_offset] = ...
114 model.implicit_gradient_operators_algorithm(model, model_data, ...
116 model.t = model.t - model.dt;
117 gradL_I = gradL_I * model.dt + speye(size(gradL_I));
118 rhs = - V(:,n) + U(:,t) - exprhs - model.dt * Urhs;
120 residual = sqrt(model.inner_product(model, model_data, VU, VU));
122 if model.newton_regularisation
123 [i] = find(Vnew < 0);
125 % epsilon = 0.99*min(-V(i,n)./VU(i));
127 disp(
' == regularisation == ');
129 %Vnew = V(:,n) + epsilon * VU;
133 nres = Vnew - U(:,t) + exprhs + model.dt * Urhs;
134 nres = sqrt(model.inner_product(model, model_data, nres, nres));
135 if model.verbose > 10
136 disp([num2str(residual),
' newton residual: ', num2str(nres)]);
138 if nres < model.newton_epsilon
144 Unewton = [Unewton, V(:,1:n+1)];
145 nlast = nlast + n + 2;
146 else % scheme with implicit parts, but without newton_scheme
147 if ~model.data_const_in_time
149 model.t=model.t+model.dt;
150 [L_diff_I, const_diff_I] = model.implicit_operators_algorithm(model, model_data, []);
151 model.t=model.t-model.dt;
154 inc = inc - const_diff_I;
156 L_I = L_I * model.dt + speye(size(L_I));
157 % L_I = speye(size(L_diff_I));
160 rhs = U(:,t) - model.dt * inc; % + model.dt * L_diff_I * U(:,t);
162 U(:,t+1) = L_I \ rhs;
163 % [U(:,t+1), flag] = bicgstab(L_I, rhs, 10e-5, 1000);
165 disp([
'error in system solution, solver return flag = ', ...
172 if model.newton_solver && model.verbose > 0
173 disp([
'computation ready: We needed at most ', num2str(nmax),
' Newton steps']);
178 if model.newton_solver
179 sim_data.Unewton = Unewton;
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 INC = fv_conv_explicit_space(U, NU_ind,gridbase grid, params)
fv_conv_explicit_space(U,NU_ind,grid,params) function applying an FV-space-discretization operator st...
function clear_all_caches()
This function clears the caches generated for caching of gradient data in evolution schemes...