rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
ldg_operators_explicit.m
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)
3 %
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
10 %
11 % required fields of model:
12 % dt : timestep
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
16 %
17 % optional fields of model:
18 % mu_names : names of fields to be regarded as parameters in vector mu
19 
20 % Bernard Haasdonk 27.8.2009
21 
22 % determine affine_decomposition_mode as integer
23 
24 disp('to be adjusted!');
25 keyboard;
26 
27 decomp_mode = model.decomp_mode;
28 
29 % call operator function with or without grid
30 if decomp_mode < 2
31  n = length(model_data.grid.A);
32 end;
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, ...
36  model_data);
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);
40 
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})))]);
49 % end
50 % keyboard;
51 %end;
52 
53 %if decomp_mode==0
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));
58 % if model.verbose>=9
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)]);
63 % end;
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)]);
67 % end;
68 % end;
69 %end;
70 
71 if decomp_mode == 2
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(:)];
80  L_I = 1;
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(:)];
90  L_I = {speye(n)};
91  L_E = [{speye(n)}; L_E_conv(:)];
92  b = bdir_E_conv(:);
93  % keyboard;
94 else % decomp_mode = 0 => complete
95  L_I = speye(n);
96  L_E = speye(n) - model.dt * L_E_conv;
97  b = model.dt * (bdir_E_conv);
98 end;
99 
100 %if decomp_mode == 0
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;
105 %
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
110 % keyboard;
111 % end;
112 % end;
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
117 % keyboard;
118 % end;
119 % end;
120 % if ~isempty(find(sum(L,2)<0))
121 % disp('warning: explicit matrix has negative row sum contributions!!');
122 % if model.verbose>=10
123 % keyboard;
124 % end;
125 % end;
126 % end;
127 %
128 % % only report this time consuming computation once:
129 % if (model.t==0) && (model.verbose >= 10)
130 % fv_estimate_operator_norms(model,model_data);
131 % end;
132 %end;
133 
134 %| \docupdate
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17