3 % Function computing a numerical convective Lax-Friedrichs
6 % Dirichlet-boundary treatment is performed, but Neumann-values are set to
7 %
'nan'. They must be performed by the calling
function.
10 % ``g_{jl}(u,v) = \frac12 (f(u,x_{jl}) + f(v,x_{jl})) \cdot \nu_{jl} -
11 % \frac{|S_{jl}|}{2 \lambda_{jl}} (v-u)``
13 % ``g_{jl}(u,v) = \frac12 |S_{jl}| \left( (f(u)+f(v)) \cdot n_{jl} -
14 % \frac{1}{\lambda_{jl}} (v-u) \right)``
15 % with `\lambda_{jl} := \frac{\Delta t}{|m_j - m_l|}`.
17 % generated fields of num_flux:
18 % Lg: Lipschitz constant satisfying
19 % ``|g_{ij}(w,v) - g_{ij}(w
',v')| \leq Lg |S_{ij}| ( |w
'-w| + |v'-v| )``
20 % G: the matrix of numerical flux values across the edges.
21 % boundary treatment is performed as follows:
22 % - in dirichlet-boundary points, the
23 % neighbour-values are simply set to
this value and
24 % the flux is computed.
25 % - Neumann-treatment is not performed here
28 % required fields of model:
29 % lxf_lambda: (inverse of) numerical diffusion in Lax-Friedrichs flux
31 grid = model_data.grid;
33 if ~isempty(model.flux_quad_degree)
34 model.flux_quad_degree = 1;
38 num_flux.G = nan* ones(size(grid.NBI));
40 % determine index sets of boundary types
41 real_NB_ind = find(grid.NBI>0);
42 NBIind = grid.NBI(real_NB_ind);
43 INBind = grid.INB(real_NB_ind);
44 % NB_real_ind comprises some kind of inverse cell-neighbour relation that
45 % leads back to the original cell from where we started:
48 % [i,j]=ind2sub(size(grid.NBI),NB_real_NB_ind(sub2ind(size(grid.NBI),k,l)))
51 NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);
53 dir_NB_ind = find(grid.NBI == -1);
55 % matrix with the neighbouring function values
56 % UNB = nan * ones(size(flux_mat.Fx));
57 UNB = nan * ones(size(grid.NBI));
58 UNB(real_NB_ind) = U(grid.NBI(real_NB_ind));
59 UU = repmat(U,1,size(grid.NBI,2));
61 % 'FxNB' and 'FyNB' and 'UNB' are computed, where possible:
62 % both are valid in 'real_NB_ind'
64 % matrix evaluating the flux by quadratures over all edges
67 % in case of divergence cleaning, the velocity field must be
68 % available!! This is obtianed by setting u to 1
69 if (model.divclean_mode>0)
70 Udummy = ones(size(U));
74 % matrix with flux values, the neighbour values U(j) inserted
75 % values must be reflected across all edges
76 % i.e. GNB(i,j) = flux.G(NBI(i,j),INB(i,j))
77 % take care of boundary-indices!!
78 FxNB = nan * ones([grid.nelements,grid.nneigh]);
79 FyNB = nan * ones([grid.nelements,grid.nneigh]);
80 % [i,j] = ind2sub([grid.nelements, grid.nneigh],NB_real_NB_ind);
81 Fx = reshape(flux_mat(1,:,:),size(FxNB));
82 Fy = reshape(flux_mat(2,:,:),size(FxNB));
83 FxNB(real_NB_ind) = Fx(NB_real_NB_ind);
84 FyNB(real_NB_ind) = Fy(NB_real_NB_ind);
86 % L1 is the supremum in space of the derivative of the first component
87 % L2 is the supremum in space of the derivative of the second component
88 L1 = max(max(abs(Fx)));
89 L2 = max(max(abs(Fy)));
91 % Lipschitz-constant of LxF Flux
92 num_flux.Lg = 1/(2*lambda) + 1/2 * norm([L1, L2]);
93 % an alternative formula without motivation:
94 % num_flux.Lg = max(max(grid.S)) * model.c * 0.25;
96 % in the linear case the following worked:
97 % numflux_mat.G = 1/2 * grid.S .* ...
98 % (grid.Nx.* (flux_mat.Vx.*(OU + OUext(grid.NBI))) + ...
99 % grid.Ny.* (flux_mat.Vy.*(OU + OUext(grid.NBI)))) ...
100 % + 1/(2 * flux_mat.lambda) * grid.S .* (OU - OUext(grid.NBI));
102 if ~isempty(model.lxf_lambda)
103 lambda = model.lxf_lambda;
106 % in the general nonlinear case:
107 num_flux.G = 1/2 * grid.EL .* ...
108 (grid.NX.* (Fx + FxNB) + ...
109 grid.NY.* (Fy + FyNB)) ...
110 + 1/(2 * lambda) * grid.EL .* (UU - UNB);
112 % Dirichlet boundary treatment
113 if ~isempty(dir_NB_ind>0)
114 % determine dirichlet-boundary values as required by convective
115 % and diffusive flux. Neumann determined at end
116 Xdir = grid.ECX(dir_NB_ind);
117 Ydir = grid.ECY(dir_NB_ind);
118 Udir = model.dirichlet_values_ptr([Xdir,Ydir],model);
120 % if conv_flux_matrix has been divergence cleaned, no further
121 % flux evaluations may be performed!!
122 if model.divclean_mode > 0
124 disp('skipping boundary flux eval. due to divergence cleaning!!');
126 % the following is unstable, if U reaches 0!!!!
127 % Fdir.Fx = flux_mat.Fx(dir_NB_ind)./UU(dir_NB_ind).*Udir;
128 % Fdir.Fy = flux_mat.Fy(dir_NB_ind)./UU(dir_NB_ind).*Udir;
129 % instead again evaluation of the whole matrix.
130 FOnex = reshape(flux_mat_UOne(1,:,:),size(FxNB));
131 FOney = reshape(flux_mat_Uone(2,:,:),size(FxNB));
132 Fdir = [FOnex(dir_NB_ind),FOney(dir_NB_ind)].*Udir;
134 else % ordinary flux evaluation by quadratures
135 % Fdir = conv_flux(model,Udir, Xdir, Ydir);
136 [elids, edgeids] = ind2sub(size(grid.VI),dir_NB_ind);
137 PP = edge_quad_points(grid,elids, edgeids, ...
138 model.flux_quad_degree);
140 % the following is only relevant in case of use of a
141 % model.use_velocitymatrix_file and filecaching mode 2
142 if isfield(model, 'velocity_matrixfile_extract') && ...
143 model.filecache_velocity_matrixfile_extract == 2;
144 model.velocity_matrixfile = ...
145 cache_velocity_matrixfile_extract(model,PP(:,1),PP(:,2),...
149 ff = model.conv_flux_ptr(PP, repmat(Udir(:),model.flux_quad_degree,1), ...
151 Fdir.Fx = edge_quad_eval_mean(grid,elids,edgeids,...
152 model.flux_quad_degree,ff(:,1));
153 Fdir.Fy = edge_quad_eval_mean(grid,elids,edgeids,...
154 model.flux_quad_degree,ff(:,2));
157 % num_flux.G(dir_NB_ind) = grid.S(dir_NB_ind) .* ...
158 % (grid.Nx(dir_NB_ind) .* Fdir.Fx + ...
159 % grid.Ny(dir_NB_ind) .* Fdir.Fy);
161 % better solution: Lxf with "ghost-cells"
162 num_flux.G(dir_NB_ind) = 1/2 * grid.EL(dir_NB_ind) .* ...
163 (grid.NX(dir_NB_ind).* (Fx(dir_NB_ind) + Fdir.Fx) + ...
164 grid.NY(dir_NB_ind).* (Fy(dir_NB_ind) + Fdir.Fy)) ...
165 + 1/(2 * lambda) * grid.EL(dir_NB_ind) .* ...
166 (UU(dir_NB_ind) - Udir);
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 [ flux_mat , lambda ] = fv_conv_flux_matrix(model, model_data, U, conv_flux_ptr)
function computing the flux matrix of a convection problem. simply reformatting the grid data suitabl...
function num_flux = fv_num_conv_flux_lax_friedrichs(model, model_data, U, NU_ind)
Function computing a numerical convective Lax-Friedrichs flux matrix.