fv_num_conv_flux_waterflow_upwind.m

function num_flux = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U)
%function num_flux = fv_num_conv_flux_engquist_osher(model,model_data,S,U)
%
% Engquist-Osher-Flux in 2D:
%   ``g_{jl}(u,v) = |e_{jl}| \cdot (c^+_{jl}(u) + c^-_{jl}(v))``
%   with
% ``c^+_{jl}(u) = c_{jl}(0) + \int_0^u   \max(c_{jl}'(s),0) ds``
% ``c^-_{jl}(u) =             \int_0^u   \min(c_{jl}'(s),0) ds``
% ``c_{jl}(u)   = f(u) \cdot n_{jl}``
%
%   if we assume, that
%         -# `f'(u)\cdot n` does not change sign on `[0,u]` and
%         -# `f(0) = 0`
%  then
%    this simplifies to
%     ``c^+_{jl}(u) = \left\{\begin{array}{rl}
%         f(u)\cdot n_{jl} & \mbox{if }f'(u)\cdot n_{jl} > 0 \\
%         0                & \mbox{otherwise}
%       \end{array}\right.``
%     ``c^-_{jl}(v) = \left\{\begin{array}{rl}
%         f(v)\cdot n_{jl} & \mbox{if }f'(v)\cdot n_{jl} < 0 \\
%         0                & \mbox{otherwise}.
%       \end{array}\right.``
%
% Dirichlet-boundary treatment is performed, but
% Neumann-values are set to 'nan'. They must be performed by
% the calling function.
%
% fields of num_flux:
%   Lg: lipschitz constant satisfying
%       ``|g_{ij}(w,v) - g_{ij}(w',v')| \leq Lg |S_{ij}| ( |w'-w| + |v'-v| ).``
%   G:  the matrix of numerical flux values across the edges.
%       boundary treatment is performed as follows:
%         - in dirichlet-boundary points, the
%         neighbour-values are simply set to this value and
%         the flux is computed.
%         - Neumann-treatment is not performed here
%         .
%
% We assume for the convective flux function:
%         -# `f'(u)\cdot n` does not change sign on interval `(u,0)`
%         -# `f(0) = 0`
%
% Additionally, a velocity function must be available in order to
% compute the derivative of the convection term.
%

% Martin Drohmann 17.10.2011

  grid = model_data.grid;

  if ~isfield(model, 'flux_quad_degree') || isempty(model.flux_quad_degree)
    model.flux_quad_degree = 1;
  end;

  num_flux = [];
  num_flux.G = nan* ones(size(grid.NBI));

  % determine index sets of boundary types
  real_NB_ind = find(grid.NBI>0);
  NBIind = grid.NBI(real_NB_ind);
  INBind = grid.INB(real_NB_ind);
  % NB_real_ind comprises comprises some kind of inverse cell-neighbour
  % relation that leads back to the original cell from where we started
  % Falls 
  %  [i,j]=ind2sub(size(grid.NBI),NB_real_NB_ind(sub2ind(size(grid.NBI),k,l)))
  % dann ist
  %  grid.NBI(i,j) == k
  NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);

  dir_NB_ind = find(grid.NBI == -1);
  neu_NB_ind = (grid.NBI == -2);

  %%%%%%%%%%%%%%%%%%%  convective numerical flux: Engquist Osher %%%%%%%%%%%%%%

  % evaluate flux matrix
  % FxNB and FyNB and UNB are computed, where possible:
  %    both are valid in real_NB_ind

  % matrix evaluating the flux by quadratures over all edges

  Sn = repmat(S, 1, size(grid.ECX,2));
  [elids, edgeids] = ind2sub(size(grid.VI),1:length(grid.VI(:)));
  PP = edge_quad_points(grid,elids, edgeids, model.flux_quad_degree);
  fs = model.two_phase.waterflow(PP,repmat(Sn(:), model.flux_quad_degree, 1)', model);

  %disp(fs(64));

  FS_NB = nan * ones([grid.nelements,grid.nneigh]);
  %  [i,j] = ind2sub([grid.nelements, grid.nneigh],NB_real_NB_ind);
  FS = reshape(fs,size(FS_NB));
  FS_NB(real_NB_ind) = FS(NB_real_NB_ind);

  UU = zeros(size(model_data.gEI));
  UU(model_data.gEI >0) = U(model_data.gEI(model_data.gEI >0));
  UU(neu_NB_ind) = 0;

  % evaluate flux matrix derivative (average over edges)
  flux_derivative_mat.Vx = grid.NX .* UU;
  flux_derivative_mat.Vy = grid.NY .* UU;
  % produced fields: Vx, Vy

  % matrix with flux values, the neighbour values U(j) inserted
  % VxNB = nan * ones(size(flux_derivative_mat.Vx));
  % VxNB(real_NB_ind) = flux_derivative_mat.Vx(NB_real_NB_ind);

  % VyNB = nan * ones(size(flux_derivative_mat.Vy));
  % VyNB(real_NB_ind) = flux_derivative_mat.Vy(NB_real_NB_ind);

  if ~isempty(dir_NB_ind>0)
    error('Dirichlet boundaries not yet implemented for waterflow updwind');
    % ordinary flux evaluation by quadratures
    Xdir = grid.ECX(dir_NB_ind);
    Ydir = grid.ECY(dir_NB_ind);
    FSdir = model.dirichlet_values_S_ptr([Xdir,Ydir],model);

    [elids, edgeids] = ind2sub(size(grid.VI),dir_NB_ind);
    PP = edge_quad_points(grid,elids,edgeids,model.flux_quad_degree);
    fs = model.two_phase.waterflow(PP,...
           repmat(FSdir(:),model.flux_quad_degree,1), ...
           model );

    Sdir = edge_quad_eval_mean(grid,elids,edgeids,...
      model.flux_quad_degree,fs(:));

    Udir = U(model_data.gEI(dir_NB_ind));
    ffdir_derivative_mat.Vx = grid.NX(dir_NB_ind) .* Udir;
    ffdir_derivative_mat.Vy = grid.NY(dir_NB_ind) .* Udir;

    Fdir_derivative.Vx = edge_quad_eval_mean(grid,elids,edgeids,...
      model.flux_quad_degree,...
      ffdir_derivative_mat.Vx);
    Fdir_derivative.Vy = edge_quad_eval_mean(grid,elids,edgeids,...
      model.flux_quad_degree,...
      ffdir_derivative_mat.Vy);

    FS_NB(dir_NB_ind)  = Sdir(:);
    % VxNB(dir_NB_ind)   = Fdir_derivative.Vx(:);
    % VyNB(dir_NB_ind)   = Fdir_derivative.Vy(:);

  end;

  % the following nicely maintains nan in Neumann-boundary edges
  c_jl_u_deri = (flux_derivative_mat.Vx + flux_derivative_mat.Vy);
  I_c_jl_u_deri_pos = find (c_jl_u_deri > 0);
  I_c_jl_v_deri_neg = c_jl_u_deri <= 0;

  G = nan * ones(size(grid.NBI)); % matrix of Nans
  G(I_c_jl_u_deri_pos) = ...
    FS(I_c_jl_u_deri_pos) .* ...
    (flux_derivative_mat.Vx(I_c_jl_u_deri_pos) + flux_derivative_mat.Vy(I_c_jl_u_deri_pos));
  G(I_c_jl_v_deri_neg) = ...
    FS_NB(I_c_jl_v_deri_neg) .* ...
    (flux_derivative_mat.Vx(I_c_jl_v_deri_neg) + flux_derivative_mat.Vy(I_c_jl_v_deri_neg));

  % add both components for final flux in all inner edges
%  num_flux.G = grid.EL.* G;
  num_flux.G = G;

  % Lipschitz konstant Lg
  %       |g_ij(w,v) - g_ij(w',v')| <= Lg |S_ij| ( |w'-w| + |v'-v| ) 

  lambda = Inf;
  num_flux.Lg = 1/ lambda;