3 % computes a neumann contribution matrix
for finite volume time evolution operators,
4 % <b> or their Frechet derivative </b>
6 % This
function computes the neumann boundary contribution
operator
7 % `L_{\text{neu}}` and a corresponding vector ` b_{\text{neu}}` that can be used
8 % by fv_operators_implicit_explicit() to build evolution matrices for a finite
9 % volume time step `L_I U^{k+1} = L_E U^k + b_E + b_I`. This operator
10 % contribution must be activated via the flag
11 % 'model.operators_neumann_implicit'
13 % The analytical terms inspiring this operator look like
14 % ` f(u) \cdot n = b_{\text{neu}}(u) \quad
15 % \text{on }\partial \Omega_{\text{neu}},`
16 % where `f` can be some flux
function, e.g. chosen via the
17 %
'model.conv_flux_ptr' and `b_{\text{neu}}(u)` is the neumann boundary
18 % condition selected via
'model.neumann_values_ptr'.
22 % fv_operators_diff_implicit_gradient_tensor()
24 % Required fields of model:
25 % neumann_values_ptr :
28 % L_E_neu : sparse matrix `L_{\text{neu}}`
29 % b_E_neu : offset vector `b_{\text{neu}}`
31 % Bernard Haasdonk 13.7.2006
33 % determine affine_decomposition_mode as integer
34 decomp_mode = model.decomp_mode;
37 if ~isempty(model_data)
38 grid = model_data.grid;
42 % only get and forward sigmas:
43 bneu0 = model.neumann_values_ptr([],[],[],model);
47 else % decomp_mode == 0 or 1
50 % determine index sets of boundary types
51 real_NB_ind = find(grid.NBI>0);
52 % dir_NB_ind =find(grid.NBI == -1);
53 neu_NB_ind =find(grid.NBI == -2);
56 disp(['found ',num2str(length(real_NB_ind)),' non-boundary edges.'])
57 disp(['found ',num2str(length(neu_NB_ind)),' neumann-boundary edges.'])
60 if ~isempty(neu_NB_ind)
61 % determine linearization of neumann boundary values
62 % FNneu = aneu * u + bneu
63 Xneu = grid.ECX(neu_NB_ind);
64 Yneu = grid.ECY(neu_NB_ind);
67 Uneu = zeros(length(Xneu),1);
68 bneu0 = model.neumann_values_ptr([Xneu(:),Yneu(:)],...
70 [grid.NX(neu_NB_ind),grid.NY(neu_NB_ind)],...
73 Uneu = ones(length(Xneu),1);
74 bneu1 = model.neumann_values_ptr([Xneu(:),Yneu(:)],...
76 [grid.NX(neu_NB_ind),grid.NY(neu_NB_ind)],...
79 %%%%%%%% explicit matrix:
81 % (L_E_neu)il = 1/|Ti| *
82 % (sum j NBneu(i) |S_ij| a_neu_ij for l=i
87 val = zeros(size(grid.ESX));
88 val(neu_NB_ind)= grid.EL(neu_NB_ind).* aneu;
89 L_E_neu = sparse(1:n,1:n, grid.Ainv(:).* sum(val,2));
92 if ~isempty(find(diag(L_E_neu)<0,1))
93 disp('warning: neumann-boundary has inflow velocity !!!! ');
97 %%%%%%%% offset-vector:
98 % (b_E_neu)_i = - 1/|T_i| * ...
99 % sum_j NBneu(i) |S_ij| b_E_neu_ij
100 val(neu_NB_ind) = grid.EL(neu_NB_ind).* bneu0;
101 b_E_neu = - grid.Ainv(:) .* sum(val,2);
102 else % decomp_mode == 1 % identical computation for all components
103 Q_bneu = length(bneu0);
104 L_E_neu = cell(Q_bneu,1);
105 b_E_neu = cell(Q_bneu,1);
107 aneu = bneu1{q}-bneu0{q};
108 val = zeros(size(grid.ESX));
109 val(neu_NB_ind)= grid.EL(neu_NB_ind).* aneu;
110 L_E_neu{q} = sparse(1:n,1:n, grid.Ainv(:).* sum(val,2));
111 val(neu_NB_ind)= grid.EL(neu_NB_ind).* bneu0{q};
112 b_E_neu{q}= - grid.Ainv(:) .* sum(val,2);
114 end; % decomp more select
115 else % compute frechet derivative ( non-linear flux )
116 old_conv_flux_ptr = model.conv_flux_ptr;
117 model.conv_flux_ptr = model.conv_flux_derivative_ptr;
119 bneu = model.neumann_values_ptr([Xneu(:),Yneu(:)], U,...
120 [grid.NX(neu_NB_ind), grid.NY(neu_NB_ind)],...
122 model.conv_flux_ptr = old_conv_flux_ptr;
124 error(
'only decomp mode == 0 is allowed!')
126 val = zeros(size(grid.ESX));
127 val(neu_NB_ind)= grid.EL(neu_NB_ind).* bneu;
128 L_E_neu = sparse(1:n,1:n, grid.Ainv(:).* sum(val,2));
131 % check nonnegativity
132 if ~isempty(find(diag(L_E_neu)<0,1))
133 disp('warning: neumann-boundary has inflow velocity !!!! ');
136 else % neumann boundary empty
138 L_E_neu = sparse(n,n);
139 b_E_neu = zeros(n,1);
140 else % decomp_mode == 1
141 % L_E_neu = {}; %{sparse(n,n)};
142 % b_E_neu = {}; % {zeros(n,1)};
143 % cannot determine in mode==2, whether neumann boundary is empty. So
144 % at least one component must be provided.
145 L_E_neu = {sparse(n,n)};
146 b_E_neu = {zeros(n,1)};
147 end; % neumann boundary nonempty
148 end; % neumann-boundary empty select
150 end; % decomp_mode ==0 or 1