rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
fv_num_conv_flux_lax_friedrichs.m
Go to the documentation of this file.
1 function num_flux = fv_num_conv_flux_lax_friedrichs(model,model_data,U,NU_ind)
2 % function num_flux = fv_num_conv_flux_lax_friedrichs(model,model_data,U,NU_ind)
3 % Function computing a numerical convective Lax-Friedrichs
4 % flux matrix.
5 %
6 % Dirichlet-boundary treatment is performed, but Neumann-values are set to
7 % 'nan'. They must be performed by the calling function.
8 %
9 % Lax-Friedrichs-Flux:
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)``
12 % or equivalently
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|}`.
16 %
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
26 % .
27 %
28 % required fields of model:
29 % lxf_lambda: (inverse of) numerical diffusion in Lax-Friedrichs flux
30 
31 grid = model_data.grid;
32 
33 if ~isempty(model.flux_quad_degree)
34  model.flux_quad_degree = 1;
35 end;
36 
37 num_flux = [];
38 num_flux.G = nan* ones(size(grid.NBI));
39 
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:
46 %
47 % If
48 % [i,j]=ind2sub(size(grid.NBI),NB_real_NB_ind(sub2ind(size(grid.NBI),k,l)))
49 % then
50 % grid.NBI(i,j) == k
51 NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);
52 
53 dir_NB_ind = find(grid.NBI == -1);
54 
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));
60 
61 % 'FxNB' and 'FyNB' and 'UNB' are computed, where possible:
62 % both are valid in 'real_NB_ind'
63 
64 % matrix evaluating the flux by quadratures over all edges
65 [flux_mat,lambda] = fv_conv_flux_matrix(model,model_data,U);
66 
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));
71  [flux_mat_Uone,lambda_one] = fv_conv_flux_matrix(model,model_data,Udummy);
72 end;
73 
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);
85 
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)));
90 
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;
95 
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));
101 
102 if ~isempty(model.lxf_lambda)
103  lambda = model.lxf_lambda;
104 end;
105 
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);
111 
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);
119 
120  % if conv_flux_matrix has been divergence cleaned, no further
121  % flux evaluations may be performed!!
122  if model.divclean_mode > 0
123  if model.verbose>9
124  disp('skipping boundary flux eval. due to divergence cleaning!!');
125  end;
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;
133 
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);
139 
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),...
146  'dirichlet_bnd');
147  end;
148 
149  ff = model.conv_flux_ptr(PP, repmat(Udir(:),model.flux_quad_degree,1), ...
150  model );
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));
155  end;
156 
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);
160 
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);
167 
168 end
169 
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
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.