fv_num_diff_flux_pressure_gradient.m

Go to the documentation of this file.

function num_flux_mat = fv_num_diff_flux_pressure_gradient(model, model_data, P, S)
%function num_flux_mat = fv_num_diff_flux_gradient(model, model_data, U)
% computes a numerical diffusive flux for a diffusion problem
%
% This function computes a numerical diffusive flux matrix for a diffusion
% problem. Dirichlet boundary treatment is performed, but neuman values are set
% to 'NaN'. These must be handled by the calling function.
%
% The analytical term inspiring this flux looks like
% `` \nabla \cdot \left(d(u) \nabla l(u) \right). ``
% Here, the functions `d(u)` and `l(u)` differ in the numerical implementation.
% The latter is evaluated in cell centers during flux computation, whereas the
% diffusion functional `d(u)` is evaluated on the cell edges by averaging the
% solutions value there.
%
% The implemented flux then can be written for inner boundaries as
%   `` g_{ij}(u,v) = d(\{u,v\}) \frac{l(v) - l(u)}{|d_{ij}|} |e_{ij}| ``
% and for Dirichlet boundaries as
%   `` g_{dir,ij}(u) = d(u_{dir}(x_{ij}))
%          \frac{l(u_{dir}(x_{ij})) - l(u)}{|d_{ij}|} |e_{ij}|, ``
% where `\{u,v\}` refers to the average of `u` and `v`.
%
% Parameters:
%   P: Dof vector of finite volume function `p`, i.e. pressure argument for
%      this operator.
%
% required fields of model:
%   diffusivity_ptr:  function pointer to diffusivity function `d`
%   laplacian_ptr:    function pointer to non-linear function `l`
%
% generated fields of num_flux_mat:
%   epsilon: diffusion coefficient bound
%   G:  the matrix of numerical flux values across the edges.
%       boundary treatment is performed as follows:
%       in dirichlet-boundary edges, the neighbour-values are simply set
%       to this value and the flux is computed.
%

% Bernard Haasdonk 10.4.2006

grid = model_data.grid;

num_flux_mat = [];
num_flux_mat.epsilon = 0;
num_flux_mat.G = zeros(model_data.gn_edges, 1);

% determine index sets of boundary types
real_NB_ind = find(grid.NBI>0 & grid.NBI < repmat((1:grid.nelements)', 1, grid.nneigh));
%  NBIind = grid.NBI(real_NB_ind);
%  INBind = grid.INB(real_NB_ind);
%  NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);

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

if model.verbose >= 29
  disp(['found ',num2str(length(real_NB_ind)),' non-boundary edges.'])
  disp(['found ',num2str(length(dir_NB_ind)),' dirichlet-boundary edges.'])
  disp(['found ',num2str(length(neu_NB_ind)),' neumann-boundary edges.'])
end;

tmpS = repmat(S, 1, grid.nneigh);
neiS = tmpS;
real_nb_ind = grid.NBI > 0;
neiS(real_nb_ind) = S(grid.NBI(real_nb_ind));
%edgeU = 0.5 * (tmpU + neiU);
%edgeU = sqrt(tmpU.*neiU);
%edgeU = min(tmpU, neiU);
%edgeU = 2*(tmpU.*neiU)./(tmpU+neiU);

difftmp = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], tmpS(:), model);
diffnei = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], neiS(:), model);

meanptr = model.mean_ptr;

diff.K       = meanptr([difftmp.K, diffnei.K]);
diff.epsilon = meanptr(max([difftmp.epsilon, diffnei.epsilon],0));
if length(diff.K) == 1
  K = ones(size(grid.ESX))*diff.K;
else
  K = reshape(diff.K,size(grid.ESX));
end
nfaces = size(grid.NBI,2);

PP = repmat(P,1,nfaces);
% matrix with the neighbouring values products
PNB = nan * ones(size(PP));
PNB(real_NB_ind) = PP(grid.NBI(real_NB_ind));

% set real neighbour values
orientation = grid.NX(real_NB_ind) + grid.NY(real_NB_ind);
num_flux_mat.G(1:model_data.gn_inner_edges) = ...
    K(real_NB_ind).*grid.EL(real_NB_ind).*(grid.DS(real_NB_ind).^(-1).* ...
                                          (PNB(real_NB_ind)- ...
                                           PP(real_NB_ind))) ...
    .* orientation;

% determine dirichlet-boundary values as required by convective
% and diffusive flux.
if ~isempty(dir_NB_ind)
  Xdir  = grid.ESX(dir_NB_ind);
  Ydir  = grid.ESY(dir_NB_ind);
  Sdir  = model.dirichlet_values_s_ptr([Xdir,Ydir],model);
  Pdir  = model.dirichlet_values_p_ptr([Xdir,Ydir],model);
  Kdir  = model.two_phase.toal_mobility([Xdir,Ydir], Sdir, model);
  Kdir  = Kdir.K;
  % determine distances circumcenter to edge-midpoint for
  % gradient approximation
  %[dir_NB_i,dir_NB_j] = ind2sub(size(UU),dir_NB_ind);
  %Ddir = sqrt((grid.ESX(dir_NB_ind)-grid.CX(dir_NB_i)).^2 + ...
    %         (grid.ESY(dir_NB_ind)-grid.CY(dir_NB_i)).^2);
  % set dirichlet neighbour values
  %num_flux_mat.G(dir_NB_ind) = ...
  %  grid.EL(dir_NB_ind).*(Ddir.^(-1).* Kdir .*(UU(dir_NB_ind)-Udir));
  orientation = grid.NX(dir_NB_ind) + grid.NY(dir_NB_ind);
  num_flux_mat.G(model_data.gn_inner_edges+1:end) =  ...
      grid.EL(dir_NB_ind).*((1.0*grid.DS(dir_NB_ind)).^(-1).* Kdir ...
                            .*(PP(dir_NB_ind)-Pdir)) ...
      .* orientation;
end

% set diffusivity bound
num_flux_mat.epsilon = diff.epsilon;

end

% vim: set sw=2 et: