fv_frechet_operators_diff_flux_pressure_gradient.m

Go to the documentation of this file.

function spm = fv_frechet_operators_diff_flux_pressure_gradient(...
    model, model_data, P, S, NU_ind)
%function spm = fv_frechet_operators_diff_flux_pressure_gradient(...
%    model, model_data, P, S, NU_ind)
% computes a jacobian of implicit non-linear convective contributions to time
% evolution matrices at a point '(P,S)'.
%
%


if nargin < 5
  NU_ind = 1:model_data.gn_edges;
end

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

real_NB_ind = 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);
%NB_real_NB_ind = NB_real_NB_ind';
[i, j]=find(real_NB_ind);
i = i';
jj = model_data.gEI(real_NB_ind);
jj = jj(:)';
dir_NB_ind = find(grid.NBI == -1);
neu_NB_ind = find(grid.NBI == -2);

spm.sn = model_data.gn_edges;
spm.sm = length(S);
[si1, subset1, dummy] = intersect(jj, NU_ind);
subset = [subset1, subset1+length(jj)];
spm.si = [si1, si1];
spm.sj = [i, NBIind'];
spm.sj = spm.sj(subset);

spm.pn = model_data.gn_edges;
spm.pm = length(P);
spm.pi = [si1, si1];
spm.pj = [i, NBIind'];
spm.pj = spm.pj(subset);

if ~isempty(dir_NB_ind)
  error('Dirichlet boundaries not yet implemented');
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);

mobi_tmp = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], tmpS(:), model);
deri_mobi_tmp = model.two_phase.total_mobility_deriv([grid.ESX(:), grid.ESY(:)], tmpS(:), model);
mobi_nei = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], neiS(:), model);
deri_mobi_nei = model.two_phase.total_mobility_deriv([grid.ESX(:), grid.ESY(:)], neiS(:), model);

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

K = meanptr(max([mobi_tmp.K, mobi_nei.K],0));
if length(K) == 1
  K = ones(size(grid.ESX))*K;
else
  K = reshape(K,size(grid.ESX));
end

d_KK = meanx(max([mobi_tmp.K, mobi_nei.K],0)).*deri_mobi_tmp.K;
if length(d_KK) == 1
  d_KK = ones(size(grid.ESX))*d_KK;
else
  d_KK = reshape(d_KK,size(grid.ESX));
end

d_KL = meany(max([mobi_tmp.K, mobi_nei.K],0)).*deri_mobi_nei.K;
if length(d_KL) == 1
  d_KL = ones(size(grid.ESX))*d_KL;
else
  d_KL = reshape(d_KL,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;

orientation_i = grid.NX(real_NB_ind) + grid.NY(real_NB_ind);
orientation_nb = grid.NX(real_NB_ind) + grid.NY(real_NB_ind);
TKL_i = grid.EL(real_NB_ind).*grid.DS(real_NB_ind).^(-1);
TKL_nb = grid.EL(real_NB_ind).*grid.DS(real_NB_ind).^(-1);


svals1 = (d_KK(real_NB_ind).* TKL_i .* (PNB(real_NB_ind) - PP(real_NB_ind)) .* orientation_i)';
svals2 = (d_KL(real_NB_ind).* TKL_nb .* (PNB(real_NB_ind) - PP(real_NB_ind)) .* orientation_nb)';

pvals1 = -(K(real_NB_ind) .* TKL_i .* orientation_i)';
pvals2 = (K(real_NB_ind) .* TKL_nb .* orientation_nb)';

spm.svals = [ svals1, svals2 ];
spm.svals = spm.svals(subset);
spm.pvals = [ pvals1, pvals2 ];
spm.pvals = spm.pvals(subset);