fv_operators_diff_implicit_gradient2.m

Go to the documentation of this file.

function spm = fv_operators_diff_implicit_gradient2(...
    model, model_data, U, NU_ind)
%function [L_I_diff, bdir_I_diff] = ...
%    fv_operators_diff_implicit_gradient2(model,model_data[, U, NU_ind])
% computes diffusion contributions to finite volume time evolution matrices,
% <b> or their Frechet derivative </b>
%
% This function computes a diffusion operator `L_{\text{diff}}` and a
% corresponding offset vector `b_{\text{diff}}` that can be used by
% fv_operators_implicit_explicit() to build evolution matrices for a finite
% volume time step
% `L_I U^{k+1} = L_E U^k + b_E + b_I`.
%
% The analytical term inspiring this operator 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.
%
% See also: fv_num_diff_flux_gradient()
%
% 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})) - u}{|d_{ij}|} |e_{ij}|, ``
% where `\{u,v\}` refers to the average of `u` and `v`.
%
% \note that this implementation only works for trivial functions `l = id`,
% because otherwise the operator would became non-linear in the space
% variable.
%
% If, however, `l(u)` is non-trivial, the flux arguments are simply
% <b>scaled</b> by `l'(u)` giving rise to the Frechet derivative
% `D L_I|_{U}` that can be used in a Newton scheme. The offset vector
% 'bdir_I_diff' can then be ignored, of course!
%
% required fields of model:
%   diffusivity_ptr          : The diffusivity function `d(u)`
%   laplacian_derivative_ptr : The derivative `l'(u)`. See above for
%                              explanation.
%   dirichlet_values_ptr     : The dirichlet function `u_{\text{dir}}`
%
% Return values:
%  L_I_diff       : a sparse matrix with diffusion contributions to `L_I`.
%  bdir_I_diff    : and offset vector containing the Dirichlet value
%                   contributions of the diffusion parts.
%

% Bernard Haasdonk 13.7.2006

% determine affine_decomposition_mode as integer
decomp_mode = model.decomp_mode;

grid = [];
if ~isempty(model_data)
  grid = model_data.grid;
end

if nargin < 4
  NU_ind = [];
end

if decomp_mode < 2
  n = grid.nelements;
  % 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_NB_ind = sub2ind(size(grid.NBI),NBIind, INBind);
  dir_NB_ind =find(grid.NBI == -1);
  real_or_dir_NB_ind = [real_NB_ind(:); dir_NB_ind(:)];

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

% currently assuming full implicit discretization of diffusivity

%%%%%%%% 1. implicit matrix:
% has entries 
% (L_I_diff)il = 1/|T_i| sum_j (NB(i) cup NB_dir(i)) 
%                             diff*|S_ij|/d_ij            if l=i
%               -diff |S_il|/d_il                         if l is NB(i)

% evaluate all edge midpoint diffusivities at time t+dt:


if decomp_mode == 2
  diff = model.diffusivity_ptr([],[],model);
else
  if nargin == 4
    tmpU = repmat(U, 1, grid.nneigh);
    neiU = tmpU;
    real_nb_ind = grid.NBI > 0;
    neiU(real_nb_ind) = U(grid.NBI(real_nb_ind));

    %edgeU = 0.5 * (tmpU + neiU);
    %edgeU = sqrt(tmpU.*neiU);
    %edgeU = min(tmpU,neiU);
    %edgeU = 2*(1./tmpU + 1./neiU).^(-1);
    %diff = model.diffusivity_ptr([grid.ESX(:), grid.ESY(:)],edgeU(:), model);
    difftmp = model.diffusivity_ptr([grid.ESX(:),grid.ESY(:)],tmpU(:),model);
    diffnei = model.diffusivity_ptr([grid.ESX(:),grid.ESY(:)],neiU(:),model);

    if ~isfield(model, 'mean_ptr')
      model.mean_ptr = @mean;
    end
    mymean = model.mean_ptr;
    %diff.K = 2*(difftmp.K + diffnei.K) ./ (difftmp.K + diffnei.K);

    diff.K = mymean([difftmp.K, diffnei.K], 2);
    diff.epsilon = mymean([difftmp.epsilon, diffnei.epsilon]);
  else
    diff = model.diffusivity_ptr([grid.ESX(:), grid.ESY(:)],[],model);
  end
  % exact diffusivity on dirichlet edges
  if (~isempty(dir_NB_ind)) & (model.decomp_mode == 0)
    Xdir = grid.ECX(dir_NB_ind);
    Ydir = grid.ECY(dir_NB_ind);
    Udir = model.dirichlet_values_ptr([Xdir(:),Ydir(:)],model);
    Kdir = model.diffusivity_ptr([Xdir(:),Ydir(:)],Udir,model);
    diff.K(dir_NB_ind) = Kdir.K(:);
  end

end;


if length(diff.K)==1
  K = ones(size(grid.ESX))*diff.K;
else
  K = reshape(diff.K,size(grid.ESX));
end
Ulaplace = model.laplacian_derivative_ptr([grid.CX(:),...
                                           grid.CY(:)],...
                                          U, model);
if length(Ulaplace) == 1
  Ulaplace = ones(size(grid.CX)) * Ulaplace;
end
Ulaplace = repmat(Ulaplace,1,grid.nneigh);
if ~isempty(dir_NB_ind)
  Udirlaplace = model.laplacian_derivative_ptr([Xdir,Ydir],...
                                              Udir, model);
else
  Udirlaplace = [];
end
% compute products
val = zeros(size(grid.ESX));
val(real_or_dir_NB_ind) = K(real_or_dir_NB_ind).* ...
    grid.EL(real_or_dir_NB_ind).* ...
    [ Ulaplace(real_NB_ind); Udirlaplace ] .* ...
    (([grid.DS(real_NB_ind);1.0*grid.DS(dir_NB_ind)]).^(-1));
% Ti_inv = repmat(grid.Ainv(:),1,grid.nneigh);
% val = val .* Ti_inv;

% diagonal values:
%L_I_diag = sparse(1:n,1:n, sum(val,2));
% off-diagonal values:
[i,dummy] = ind2sub(size(grid.NBI),real_NB_ind);
i = i';
%L_I_offdiag = sparse(i,grid.NBI(real_NB_ind), -val(real_NB_ind),n,n);
spm.sn = n;
spm.sm = n;
spm.si = [1:n, i];
spm.sj = [1:n, grid.NBI(real_NB_ind)'];
spm.svals = [sum(val,2)', -val(real_NB_ind(:)')];

if ~isempty(dir_NB_ind)
  error('Dirichlet values not yet supported');
  % evaluate dirichlet values at edge midpoints at t+dt
  % Xdir, Ydir, Udir, K(dir_NB_ind) already computed above

  Kdir = K(dir_NB_ind); % diffusivity already computed above
  val2 = zeros(size(grid.ESX));
  val2(dir_NB_ind)= - grid.EL(dir_NB_ind).*...
    Udirlaplace.*((1.0*grid.DS(dir_NB_ind)).^(-1)) ...
    .* Kdir .*Udir;
%  bdir_I_diff = - grid.Ainv(:).* sum(val2,2);
  bdir_I_diff = - sum(val2,2);
else
  bdir_I_diff = zeros(grid.nelements,1);
end