1 function num_flux = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U)
4 % Engquist-Osher-Flux in 2D:
5 % ``g_{jl}(u,v) = |e_{jl}| \cdot (c^+_{jl}(u) + c^-_{jl}(v))``
7 % ``c^+_{jl}(u) = c_{jl}(0) + \int_0^u \max(c_{jl}
'(s),0) ds``
8 % ``c^-_{jl}(u) = \int_0^u \min(c_{jl}'(s),0) ds``
9 % ``c_{jl}(u) = f(u) \cdot n_{jl}``
12 % -# `f'(u)\cdot n` does not change sign on `[0,u]` and
16 % ``c^+_{jl}(u) = \left\{\begin{array}{rl}
17 % f(u)\cdot n_{jl} & \mbox{
if }f
'(u)\cdot n_{jl} > 0 \\
18 % 0 & \mbox{otherwise}
19 % \end{array}\right.``
20 % ``c^-_{jl}(v) = \left\{\begin{array}{rl}
21 % f(v)\cdot n_{jl} & \mbox{if }f'(v)\cdot n_{jl} < 0 \\
22 % 0 & \mbox{otherwise}.
23 % \end{array}\right.``
25 % Dirichlet-boundary treatment is performed, but
26 % Neumann-values are set to 'nan'. They must be performed by
27 % the calling function.
30 % Lg: lipschitz constant satisfying
31 % ``|g_{ij}(w,v) - g_{ij}(w
',v')| \leq Lg |S_{ij}| ( |w
'-w| + |v'-v| ).``
32 % G: the matrix of numerical flux values across the edges.
33 % boundary treatment is performed as follows:
34 % - in dirichlet-boundary points, the
35 % neighbour-values are simply set to this value and
36 % the flux is computed.
37 % - Neumann-treatment is not performed here
40 % We assume for the convective flux function:
41 % -# `f'(u)\cdot n` does not change sign on interval `(u,0)`
44 % Additionally, a velocity function must be available in order to
45 % compute the derivative of the convection term.
48 % Martin Drohmann 17.10.2011
50 grid = model_data.grid;
52 if ~isfield(model,
'flux_quad_degree') || isempty(model.flux_quad_degree)
53 model.flux_quad_degree = 1;
57 num_flux.G = nan* ones(size(grid.NBI));
59 % determine index sets of boundary types
60 real_NB_ind = find(grid.NBI>0);
61 NBIind = grid.NBI(real_NB_ind);
62 INBind = grid.INB(real_NB_ind);
63 % NB_real_ind comprises comprises some kind of inverse cell-neighbour
64 % relation that leads back to the original cell from where we started
66 % [i,j]=ind2sub(size(grid.NBI),NB_real_NB_ind(sub2ind(size(grid.NBI),k,l)))
69 NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);
71 dir_NB_ind = find(grid.NBI == -1);
72 neu_NB_ind = (grid.NBI == -2);
74 %%%%%%%%%%%%%%%%%%% convective numerical flux: Engquist Osher %%%%%%%%%%%%%%
76 % evaluate flux matrix
77 % FxNB and FyNB and UNB are computed, where possible:
78 % both are valid in real_NB_ind
80 % matrix evaluating the flux by quadratures over all edges
82 Sn = repmat(S, 1, size(grid.ECX,2));
83 [elids, edgeids] = ind2sub(size(grid.VI),1:length(grid.VI(:)));
84 PP = edge_quad_points(grid,elids, edgeids, model.flux_quad_degree);
85 fs = model.two_phase.waterflow(PP,repmat(Sn(:), model.flux_quad_degree, 1)', model);
89 FS_NB = nan * ones([grid.nelements,grid.nneigh]);
90 % [i,j] = ind2sub([grid.nelements, grid.nneigh],NB_real_NB_ind);
91 FS = reshape(fs,size(FS_NB));
92 FS_NB(real_NB_ind) = FS(NB_real_NB_ind);
94 UU = zeros(size(model_data.gEI));
95 UU(model_data.gEI >0) = U(model_data.gEI(model_data.gEI >0));
98 % evaluate flux matrix derivative (average over edges)
99 flux_derivative_mat.Vx = grid.NX .* UU;
100 flux_derivative_mat.Vy = grid.NY .* UU;
101 % produced fields: Vx, Vy
103 % matrix with flux values, the neighbour values U(j) inserted
104 % VxNB = nan * ones(size(flux_derivative_mat.Vx));
105 % VxNB(real_NB_ind) = flux_derivative_mat.Vx(NB_real_NB_ind);
107 % VyNB = nan * ones(size(flux_derivative_mat.Vy));
108 % VyNB(real_NB_ind) = flux_derivative_mat.Vy(NB_real_NB_ind);
110 if ~isempty(dir_NB_ind>0)
111 error('Dirichlet boundaries not yet implemented for waterflow updwind');
112 % ordinary flux evaluation by quadratures
113 Xdir = grid.ECX(dir_NB_ind);
114 Ydir = grid.ECY(dir_NB_ind);
115 FSdir = model.dirichlet_values_S_ptr([Xdir,Ydir],model);
117 [elids, edgeids] = ind2sub(size(grid.VI),dir_NB_ind);
118 PP = edge_quad_points(grid,elids,edgeids,model.flux_quad_degree);
119 fs = model.two_phase.waterflow(PP,...
120 repmat(FSdir(:),model.flux_quad_degree,1), ...
123 Sdir = edge_quad_eval_mean(grid,elids,edgeids,...
124 model.flux_quad_degree,fs(:));
126 Udir = U(model_data.gEI(dir_NB_ind));
127 ffdir_derivative_mat.Vx = grid.NX(dir_NB_ind) .* Udir;
128 ffdir_derivative_mat.Vy = grid.NY(dir_NB_ind) .* Udir;
130 Fdir_derivative.Vx = edge_quad_eval_mean(grid,elids,edgeids,...
131 model.flux_quad_degree,...
132 ffdir_derivative_mat.Vx);
133 Fdir_derivative.Vy = edge_quad_eval_mean(grid,elids,edgeids,...
134 model.flux_quad_degree,...
135 ffdir_derivative_mat.Vy);
137 FS_NB(dir_NB_ind) = Sdir(:);
138 % VxNB(dir_NB_ind) = Fdir_derivative.Vx(:);
139 % VyNB(dir_NB_ind) = Fdir_derivative.Vy(:);
143 % the following nicely maintains nan in Neumann-boundary edges
144 c_jl_u_deri = (flux_derivative_mat.Vx + flux_derivative_mat.Vy);
145 I_c_jl_u_deri_pos = find (c_jl_u_deri > 0);
146 I_c_jl_v_deri_neg = c_jl_u_deri <= 0;
148 G = nan * ones(size(grid.NBI)); % matrix of Nans
149 G(I_c_jl_u_deri_pos) = ...
150 FS(I_c_jl_u_deri_pos) .* ...
151 (flux_derivative_mat.Vx(I_c_jl_u_deri_pos) + flux_derivative_mat.Vy(I_c_jl_u_deri_pos));
152 G(I_c_jl_v_deri_neg) = ...
153 FS_NB(I_c_jl_v_deri_neg) .* ...
154 (flux_derivative_mat.Vx(I_c_jl_v_deri_neg) + flux_derivative_mat.Vy(I_c_jl_v_deri_neg));
156 % add both components for final flux in all inner edges
157 % num_flux.G = grid.EL.* G;
160 % Lipschitz konstant Lg
161 % |g_ij(w,v) - g_ij(w',v')| <= Lg |S_ij| ( |w'-w| + |v'-v| )
164 num_flux.Lg = 1/ lambda;