1 function [L_I,L_E,b] = ldg_operators_explicit(model,model_data)
2 %
function [L_I,L_E,b] = ldg_operators_explicit(model,model_data)
4 %
function computing the time evolution matrices
for a ldg time
5 % step L_I * Unew = L_E * U + b
6 % result are two sparse matrices L_I, L_E and an offset vector b
7 % Function supports affine decomposition, i.e. different operation modes
8 % guided by optional field decomp_mode in model. See also the
9 % contents.txt
for general explanation
11 % required fields of model:
13 % t : current time, i.e. discretization is taking place in t and
14 % t+dt
for explicit and implicit parts respectively
15 %
verbose : integer indicating the verbosity level
17 % optional fields of model:
18 % mu_names : names of fields to be regarded as parameters in vector mu
20 % Bernard Haasdonk 27.8.2009
22 % determine affine_decomposition_mode as integer
24 disp(
'to be adjusted!');
27 decomp_mode = model.decomp_mode;
29 % call
operator function with or without grid
31 n = length(model_data.grid.A);
33 % [L_I_diff, bdir_I_diff] = fv_operators_diff_implicit(model,model_data);
34 % [L_E_diff, bdir_E_diff] = fv_operators_diff_explicit(model,model_data);
35 [L_E_conv, bdir_E_conv] = ldg_operators_adv_explicit(model, ...
37 % [L_I_conv, bdir_I_conv] = ldg_operators_conv_implicit(model,model_data);
38 % [L_E_neu, bneu_E] = fv_operators_neuman_explicit(model,model_data);
39 % [L_I_neu, bneu_I] = fv_operators_neuman_implicit(model,model_data);
41 % check that offline and online simulation generate identical
42 % number of components,
for this run routine twice and compare numbers:
43 %
if model.verbose>=10 && decomp_mode >=1
44 % vnames = {
'L_I_diff',
'bdir_I_diff',
'L_E_diff',
'bdir_E_diff',...
45 %
'L_E_conv',
'bdir_E_conv',
'L_I_conv',
'bdir_I_conv',...
46 %
'L_E_neu',
'bneu_E',
'L_E_neu',
'bneu_I'};
47 %
for v = 1:length(vnames)
48 % disp([vnames{v},
': ',num2str(size(eval(vnames{v})))]);
54 % % check
for condition: 1-dt * L_E_conv(i,i) >= 0 as
this is the
55 % % coefficient of u_i^k in the convex-combination representation of u_i^(k+1)
56 % % maximum possible dt_max = min ( 1/L_E_conv(i,i))
57 % dt_max = min(full(diag(L_E_conv)).^(-1));
59 %
if model.dt > dt_max
60 % disp([
'current dt = ',num2str(model.dt),...
61 %
' is larger than cfl-condition!']);
62 % disp([
'cfl: dt_max = ',num2str(dt_max)]);
64 %
if model.dt < dt_max * 0.1;
65 % disp([
'current dt = ',num2str(model.dt),
' can be chosen much larger!']);
66 % disp([
'cfl: dt_max = ',num2str(dt_max)]);
72 % decomp_mode == 2 -> sigmas are required
73 % L_I = [1; model.dt* L_I_diff(:); ...
74 % model.dt * L_I_conv(:); model.dt * L_I_neu(:)];
75 % L_E = [1; - model.dt* L_E_diff(:); ...
76 % -model.dt * L_E_conv(:); -model.dt * L_E_neu(:)];
77 % b = model.dt * [ bdir_E_diff(:) ; bdir_I_diff(:); ...
78 % bdir_E_conv(:) ; bdir_I_conv(:); ...
79 % bneu_E(:); bneu_I(:)];
81 L_E = [1; -model.dt * L_E_conv(:)];
82 b = model.dt * bdir_E_conv(:);
83 elseif decomp_mode == 1
84 % in
case of components: concatenate eye and the other lists
85 % L_I = [{speye(n)}; L_I_diff(:); L_I_conv(:); L_I_neu(:)];
86 % L_E = [{speye(n)}; L_E_diff(:); L_E_conv(:); L_E_neu(:)];
87 % b = [bdir_E_diff(:); bdir_I_diff(:); ...
88 % bdir_E_conv(:); bdir_I_conv(:); ...
89 % bneu_E(:) ; bneu_I(:)];
91 L_E = [{speye(n)}; L_E_conv(:)];
94 else % decomp_mode = 0 => complete
96 L_E = speye(n) - model.dt * L_E_conv;
97 b = model.dt * (bdir_E_conv);
101 % % check L_E contributions:
102 % % matrix L should have positive diagonal and nonpositive off-diagonal
103 % % with row-sum between 0 and 1
104 % L = L_E_diff + L_E_conv + L_E_neu;
106 %
if model.verbose>=10
107 %
if ~isempty(find(diag(L)<0))
108 % disp(
'warning: explicit matrix has negative diagonal entries!!');
109 %
if model.verbose>=10
113 % Lo = L-diag(diag(L));
114 %
if ~isempty(find(Lo>0))
115 % disp(
'warning: explicit matrix has positive off-diagona entries!!');
116 %
if model.verbose>=10
120 %
if ~isempty(find(sum(L,2)<0))
121 % disp(
'warning: explicit matrix has negative row sum contributions!!');
122 %
if model.verbose>=10
128 % % only report
this time consuming computation once:
129 %
if (model.t==0) && (model.verbose >= 10)
130 % fv_estimate_operator_norms(model,model_data);
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...