fv_operators_conv_flux_waterflow_upwind.m

Go to the documentation of this file.

function [L_E_conv, bdir_E_conv] = fv_frechet_operators_conv_flux_waterflow_upwind(...
    model,model_data, S, U, S_ind, U_ind)
%function [L_E_conv,bdir_E_conv] = fv_frechet_operators_conv_flux_waterflow_upwind(...
%    model,model_data, P, S, P_ind, S_ind)
% computes convection contribution to finite volume time evolution matrices,
% <b> or their Frechet derivative </b>
%


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

real_NB_ind = find(grid.NBI>0);
NBIind = grid.NBI(real_NB_ind);
INBind = grid.INB(real_NB_ind);
% NB_real_ind comprises comprises some kind of inverse cell-neighbour
% relation that leads back to the original cell from where we started
% Falls 
%  [i,j]=ind2sub(size(grid.NBI),NB_real_NB_ind(sub2ind(size(grid.NBI),k,l)))
% dann ist
%  grid.NBI(i,j) == k
NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);

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


spm.sn = length(S);
spm.si;
spm.sj;
spm.snz;
spm.svals;

spm.un = length(U);
spm.ui;
spm.uj;
spm.unz;
spm.uvals;


% matrix evaluating the flux by quadratures over all edges
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);
dfs = model.two_phase.waterflow_derivative(PP,repmat(Sn(:), model.flux_quad_degree, 1)', model);
FS = reshape(fs,size(grid.NBI));
DFS = reshape(dfs, size(grid.NBI));


spm.snz =