fv_frechet_operators_conv_flux_waterflow_upwind.m

Go to the documentation of this file.

function spm = fv_frechet_operators_conv_flux_waterflow_upwind(...
    model, model_data, S, U, NU_ind)
% function spm = fv_frechet_operators_conv_flux_waterflow_upwind(...
%    model, model_data, S, U [,S_ind, U_ind])
% computes a jacobian of implicit non-linear convective contributions to time
% evolution matrices at a point '(S,U)'.
%
%


spm.sn = length(S);
spm.sm = length(S);

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

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

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

NBI = grid.NBI;
NU_ind_inv = setdiff(1:spm.sn, NU_ind);
NBI(NU_ind_inv, :) = 0;
INB = grid.INB;

real_NB_ind = NBI>0;
dir_NB_ind = find(NBI == -1);
neu_NB_ind = (NBI == -2);
%real_dir_NB_ind = [find(grid.NBI > repmat((1:grid.nelements)',1,4)); find(grid.NBI < 0)];
NBIind = NBI(real_NB_ind);
INBind = INB(real_NB_ind);
NB_real_NB_ind = sub2ind(size(NBI),NBIind,INBind);

if ~isempty(dir_NB_ind)
  error('Dirichlet boundaries not yet implemented');
end

[i,j]=find(NBI>0);
i = i';
jj = grid.NBI(real_NB_ind);
jj = jj(:)';
% Note: find(real_NB_ind) == sub2ind(size(grid.NBI),i,j));


spm.si = [NU_ind, i];
spm.sj = [NU_ind, jj];

spm.un = length(S);
spm.um = length(U);
spm.ui = i;
spm.uj = model_data.gEI(real_NB_ind);

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);
% matrix evaluating the flux by quadratures over all edges
FS  = reshape(fs,size(grid.NBI));
FS_NB = nan * ones([grid.nelements,grid.nneigh]);
FS_NB(real_NB_ind) = FS(NB_real_NB_ind);

dfs = model.two_phase.waterflow_derivative(PP,repmat(Sn(:), model.flux_quad_degree, 1)', model);
% matrix evaluating the flux derivative by quadratures over all edges
DFS = reshape(dfs,size(grid.NBI));
DFS_NB = nan * ones([grid.nelements,grid.nneigh]);
DFS_NB(real_NB_ind) = DFS(NB_real_NB_ind);

UU = nan(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

c_jl_u_deri       = (flux_derivative_mat.Vx + flux_derivative_mat.Vy);
I_c_jl_u_deri_pos = c_jl_u_deri > 0;
I_c_jl_v_deri_neg = c_jl_u_deri <= 0;
G1 = zeros(size(grid.NBI)); % zero matrix
G1(I_c_jl_u_deri_pos) = ... %grid.EL(I_c_jl_u_deri_pos) .* ...
  DFS(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));

svals1 = sum(G1, 2);

G2 = zeros(size(grid.NBI)); % zero matrix
G2(I_c_jl_v_deri_neg) = ... % grid.EL(I_c_jl_v_deri_neg) .* ...
  DFS_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));

svals2 = G2(real_NB_ind);

spm.svals = [svals1(NU_ind)', svals2'];

G3 = zeros(size(grid.NBI)); % zero matrix
G3(I_c_jl_u_deri_pos) = FS(I_c_jl_u_deri_pos) .* (grid.NX(I_c_jl_u_deri_pos) + grid.NY(I_c_jl_u_deri_pos));
G3(I_c_jl_v_deri_neg) = FS_NB(I_c_jl_v_deri_neg) .* (grid.NX(I_c_jl_v_deri_neg) + grid.NY(I_c_jl_v_deri_neg));

%G3 = grid.EL .* G3;

spm.uvals = G3(real_NB_ind);