fv_frechet_operators_diff_implicit_gradient3.m

Go to the documentation of this file.

function spm = fv_frechet_operators_diff_implicit_gradient3(...
    model, model_data, U, NU_ind)
%function spm = fv_frechet_operators_diff_implicit_gradient3(...
%    model, model_data, U)
% computes a jacobian of implicit non-linear diffusion contributions to time
% evolution matrices at a point 'U'.
%
% function computing the jacobian in 'U' of the implicit diffusion contribution
% of `L_I` and `b_I` to the time evolution matrices for a finite volume time
% step
% `L_I U^{k+1} = L_E U^k + b_E + b_I`.
% With the help of the returned Jacobian 'L_I_diff_jac' the Frechet derivative
% `DL_I ({U})` is approximated.
%
% \note The *_implicit functions perform a 'dt' increase in 'model' \em before
% evaluating the data functions.
%
% Return values:
%  spm : a struct containing all arguments for the creation of a sparse matrix
%        with the 'sparse' command.
%

% determine affine_decomposition_mode as integer

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

meanptr = model.mean_ptr;

meanx = model.mean_deriv1_ptr;
meany = model.mean_deriv2_ptr;

n = grid.nelements;

spm.sn = n;
spm.sm = n;
if nargin < 4 || isempty(NU_ind)
  NU_ind = (1:n);
end
if size(NU_ind,2)==1
  NU_ind = NU_ind';
end

NBI = grid.NBI;
NU_ind_inv = setdiff(1:n, NU_ind);
NBI(NU_ind_inv, :) = 0;
% determine index sets of boundary types
real_NB_ind = find(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(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

tmpU = repmat(U, 1, grid.nneigh);
neiU = tmpU;
real_nb_ind = NBI > 0;
neiU(real_nb_ind) = U(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);
deri_difftmp = model.diffusivity_derivative_ptr([grid.ESX(:),grid.ESY(:)],tmpU(:),model);
diffnei = model.diffusivity_ptr([grid.ESX(:),grid.ESY(:)],neiU(:),model);
deri_diffnei = model.diffusivity_derivative_ptr([grid.ESX(:),grid.ESY(:)],neiU(:),model);

meanarg = [difftmp.K, diffnei.K];
KK   = meanptr(meanarg);

nulls = sum(abs(meanarg),2)==0;
d_KK = meanx(meanarg).*deri_difftmp.K;
d_KK(nulls) = 0;
d_KL = meany(meanarg).*deri_diffnei.K;
d_KL(nulls) = 0;

diff.epsilon = meanptr([difftmp.epsilon, diffnei.epsilon]);

% exact diffusivity on dirichlet edges
if (~isempty(dir_NB_ind))
  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);
  K_deriv_dir = model.diffusivity_derivative_ptr([Xdir(:),Ydir(:)],Udir,model);
  KK(dir_NB_ind)   = Kdir.K(:);
  d_KL(dir_NB_ind) = K_deriv_dir.K(:);
  d_KK(dir_NB_ind) = 0;
end


if (length(KK) == 1)
  KK = ones(size(grid.ESX))*KK;
  d_KK = ones(size(grid.ESX))*d_KK;
  d_KL = ones(size(grid.ESX))*d_KL;
else
  KK = reshape(KK,size(grid.ESX));
  d_KK = reshape(d_KK,size(grid.ESX));
  d_KL = reshape(d_KL,size(grid.ESX));
end

Ulaplace = model.laplacian_ptr([grid.CX(:),...
                                grid.CY(:)],...
                                U, model);
DUlaplace = model.laplacian_derivative_ptr([grid.CX(:),...
                                            grid.CY(:)],...
                                           U, model);
if length(DUlaplace) == 1
  DUlaplace = ones(size(grid.CX)) * DUlaplace;
end
Ulaplace_NB = zeros(size(grid.ESX));
Ulaplace_NB(real_NB_ind) = Ulaplace(NBI(real_NB_ind));
DUlaplace_NB = zeros(size(grid.ESX));
DUlaplace_NB(real_NB_ind) = DUlaplace(NBI(real_NB_ind));

if ~isempty(dir_NB_ind)
  Ulaplace_NB(dir_NB_ind)  = model.laplacian_ptr([Xdir,Ydir],...
                                     Udir, model);
  DUlaplace_NB(dir_NB_ind) = model.laplacian_derivative_ptr([Xdir,Ydir],...
                                                Udir, model);
end

% compute products
valdiag = zeros(size(grid.ESX));

[i,dummy] = ind2sub(size(NBI),real_or_dir_NB_ind);
i = i';

valdiag(real_or_dir_NB_ind) = ...
  grid.EL(real_or_dir_NB_ind).* ...
  grid.DS(real_or_dir_NB_ind).^(-1).* ...
  ( d_KK(real_or_dir_NB_ind).* ...
      (Ulaplace(i) - Ulaplace_NB(real_or_dir_NB_ind)) ...
    + KK(real_or_dir_NB_ind).* DUlaplace(i) ...
  );

valdiag = sum(valdiag, 2);

valneigh = ...
  grid.EL(real_or_dir_NB_ind).* ...
  grid.DS(real_or_dir_NB_ind).^(-1).* ...
  ( d_KL(real_or_dir_NB_ind).* ...
      (Ulaplace(i) - Ulaplace_NB(real_or_dir_NB_ind)) ...
    - KK(real_or_dir_NB_ind).* DUlaplace_NB(real_or_dir_NB_ind) ...
  );

% 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:

spm.si = [NU_ind, i];
spm.sj = [NU_ind, NBI(real_or_dir_NB_ind)'];
spm.svals = [valdiag(NU_ind)', valneigh'];

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