rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
fv_operators_implicit_explicit_without_timestep.m
1 function [L_I,L_E,b] = fv_operators_implicit_explicit_without_timestep(model,model_data, NU_ind)
2 %function [L_I,L_E,b] = fv_operators_implicit_explicit_without_timestep(model,model_data[, NU_ind])
3 %
4 % function computing the time evolution matrices for a finite volume
5 % time step @f$ \frac{Unew-U}{\Delta t} = L_E * U + b
6 % result are two sparse matrices 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 % verbose : integer indicating the verbosity level
13 %
14 % optional fields of model:
15 % mu_names : names of fields to be regarded as parameters in vector mu
16 % decomp_mode: operation mode of the function
17 % - 0='complete' (default): no parameter dependence or decomposition is
18 % performed. output is as described above.
19 % - 1='components': For each output argument a cell array of output
20 % arguments is returned representing the q-th component
21 % independent of the parameters given in mu_names
22 % - 2='coefficients': For each output argument a cell array of output
23 % arguments is returned representing the q-th coefficient
24 % dependent of the parameters given in mu_names
25 %
26 % In 'coefficient' mode, the grid is empty
27 
28 % Markus Dihlmann 04.06.2010
29 
30 % determine affine_decomposition_mode as integer
31 decomp_mode = model.decomp_mode;
32 
33 grid = [];
34 if ~isempty(model_data)
35  grid = model_data.grid;
36 end
37 
38 if nargin < 3
39  NU_ind = [];
40 end
41 
42 % call operator function with or without grid
43 if decomp_mode < 2
44  n = length(grid.A);
45  [L_I_diff, bdir_I_diff] = model.operators_diff_implicit(model,model_data, NU_ind);
46  [L_E_diff, bdir_E_diff] = model.operators_diff_explicit(model,model_data);
47  [L_E_conv, bdir_E_conv] = model.operators_conv_explicit(model,model_data);
48  [L_I_conv, bdir_I_conv] = model.operators_conv_implicit(model,model_data);
49  [L_E_neu, bneu_E] = model.operators_neumann_explicit(model,model_data);
50  [L_I_neu, bneu_I] = model.operators_neumann_implicit(model,model_data);
51 else
52  [L_E_diff, bdir_E_diff] = model.operators_diff_explicit(model,model_data);
53  [L_I_diff, bdir_I_diff] = model.operators_diff_implicit(model,model_data, NU_ind);
54  [L_E_conv, bdir_E_conv] = model.operators_conv_explicit(model,model_data);
55  [L_I_conv, bdir_I_conv] = model.operators_conv_implicit(model,model_data);
56  [L_E_neu, bneu_E] = model.operators_neumann_explicit(model,model_data);
57  [L_I_neu, bneu_I] = model.operators_neumann_implicit(model,model_data);
58 end;
59 
60 if model.debug && decomp_mode == 0
61  % compute one assembly of affine param decomp with complete mode:
62  disp('test of affine parameter dependence of operators:')
63  success = test_affine_decomp(model.operators_diff_implicit,...
64  2,1,model,model_data);
65  success = test_affine_decomp(model.operators_diff_explicit,...
66  2,1,model,model_data);
67  success = test_affine_decomp(model.operators_conv_explicit,...
68  2,1,model,model_data);
69  success = test_affine_decomp(model.operators_conv_implicit,...
70  2,1,model,model_data);
71  success = test_affine_decomp(model.operators_neumann_implicit,...
72  2,1,model,model_data);
73  success = test_affine_decomp(model.operators_neumann_explicit,...
74  2,1,model,model_data);
75 end;
76 
77 % check that offline and online simulation generate identical
78 % number of components, for this run routine twice and compare numbers:
79 if model.verbose>=10 && decomp_mode >=1
80  disp('in fv_operators_implicit_explicit:')
81  disp(['decomp_mode = ',model.decomp_mode,', sizes of components:'])
82  vnames = {'L_I_diff','bdir_I_diff','L_E_diff','bdir_E_diff',...
83  'L_E_conv','bdir_E_conv','L_I_conv','bdir_I_conv',...
84  'L_E_neu','bneu_E','L_E_neu','bneu_I'};
85  for v = 1:length(vnames)
86  disp([vnames{v},': ',num2str(size(eval(vnames{v})))]);
87  end
88  keyboard;
89 end;
90 
91 %No implicit operators-->
92 %[L_I,b_dummy] = fv_operators_zero(model,model_data, NU_ind);
93 
94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
95 if decomp_mode == 0
96  L_I = (L_I_diff + L_I_conv + L_I_neu);
97  if(isempty(L_I))
98  L_I=sparse(n,n);
99  end
100  L_E = -(L_E_diff + L_E_conv + L_E_neu);
101  % b = model.dt * (bdir_diff + bdir_conv + bneu);
102  b = (bdir_I_diff + bdir_E_diff + bdir_I_conv + bdir_E_conv + ...
103  bneu_I + bneu_E);
104 elseif decomp_mode == 1
105  % in case of components: concatenate eye and the other lists
106  L_I = [L_I_diff(:); L_I_conv(:); L_I_neu(:)];
107  %debug, 28.10.10 --> auskommentiert
108  %if(isempty(L_I))
109  %L_I=cell(1);
110  % L_I{1} = sparse(n,n);
111  %end
112  %end debug 28.10.10
113  L_E = [L_E_diff(:); L_E_conv(:); L_E_neu(:)];
114  b = [bdir_E_diff(:); bdir_I_diff(:); ...
115  bdir_E_conv(:); bdir_I_conv(:); ...
116  bneu_E(:) ; bneu_I(:)];
117 % keyboard;
118 else
119  % decomp_mode == 2 -> sigmas are required
120  L_I = [-L_I_diff(:); -L_I_conv(:); -L_I_neu(:)];
121 %debug, 28.10.10 --> auskommentiert
122 % if(isempty(L_I))
123 % L_I=0;
124 % end
125 %end debug 28.10.10
126  L_E = [-L_E_diff(:); -L_E_conv(:); -L_E_neu(:)];
127  b = [ bdir_E_diff(:) ; bdir_I_diff(:); bdir_E_conv(:) ; bdir_I_conv(:); ...
128  bneu_E(:); bneu_I(:)];
129 % if model.t>2
130 % bdir_E_conv
131 % keyboard;
132 % end;
133 end;
134 %keyboard;
135 
136 if decomp_mode == 0
137  % check L_E contributions:
138  % matrix L should have positive diagonal and nonpositive off-diagonal
139  % with row-sum between 0 and 1
140  L = L_E_diff + L_E_conv + L_E_neu;
141 
142  if model.verbose>=10
143  if ~isempty(find(diag(L)<0))
144  disp('warning: explicit matrix has negative diagonal entries!!');
145  if model.verbose>=10
146  keyboard;
147  end;
148  end;
149  Lo = L-diag(diag(L));
150  if ~isempty(find(Lo>0))
151  disp('warning: explicit matrix has positive off-diagona entries!!');
152  if model.verbose>=10
153  keyboard;
154  end;
155  end;
156  if ~isempty(find(sum(L,2)<0))
157  disp('warning: explicit matrix has negative row sum contributions!!');
158  if model.verbose>=10
159  keyboard;
160  end;
161  end;
162  end;
163 
164  % only report this time consuming computation once:
165  if (model.t==0) && (model.verbose >= 10)
166  fv_estimate_operator_norms(model,model_data);
167  end;
168 end;
169 
170 %| \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