fv_two_phase_space.m

function [INC, J] = fv_two_phase_space(model, model_data, CU, NU_ind)
%function INC = fv_two_phase_space(model, model_data, CU, NU_ind)

  grid = model_data.grid;

  if nargin <= 3 || isempty(NU_ind)
    NU_ind = cell(1,3);
    NU_ind{1} = (1:grid.nelements)';
    NU_ind{2} = (1:model_data.gn_edges)';
    NU_ind{3} = (1:grid.nelements)';
  end

  Slength = grid.nelements;
  Ulength = model_data.gn_edges;
  S = CU(1:Slength);
  U = CU(Slength+1:Slength+Ulength);
  P = CU(Slength+Ulength+1:end);

%  model.diffusivity_ptr            = @(glob, U, model) struct('K', 1, 'epsilon', 0);%model.two_phase.phi(glob, U, model);
%  model.diffusivity_derivative_ptr = @(glob, U, model) struct('K', 0, 'epsilon', 0);%model.two_phase.phi_derivative(glob, U, model);
%  model.laplacian_ptr              = @(glob, U, model) -getfield(model.two_phase.Phi(glob, U, model), 'K')';
%  model.laplacian_derivative_ptr   = @(glob, U, model) -getfield(model.two_phase.phi(glob, U, model), 'K');

  model.diffusivity_ptr            = @(glob, U, model) model.two_phase.phi(glob, U, model);
  model.diffusivity_derivative_ptr = @(glob, U, model) model.two_phase.phi_derivative(glob, U, model);
  model.laplacian_ptr              = @(glob, U, descr) U;
  model.laplacian_derivative_ptr   = @(glob, U, descr) ones(length(U),1);

  num_diff_flux_s = fv_num_diff_flux_gradient(model, model_data, S, NU_ind{1});

  num_conv_flux_s = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind{1}, NU_ind{2});

  % homogeneous neumann boundaries...
  num_conv_flux_s.G(isnan(num_conv_flux_s.G)) = 0;

  num_diff_flux_u = fv_num_diff_flux_pressure_gradient(model, model_data, P, S, NU_ind{3}, NU_ind{1});

  tmppar.local_mass_matrix  = @fv_local_mass_matrix_rect;
  tmppar.element_quadrature = @cubequadrature;
  tmppar.evaluate_basis     = @fv_evaluate_basis;
  tmppar.pdeg = 0;
  S_upper = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.upper_source(elemin, loc, grid), 2, model_data.grid, tmppar);
  S_lower = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.lower_source(elemin, loc, grid), 2, model_data.grid, tmppar);

%  S_lower = grid.A(NU_ind{1}) .* model.two_phase.lower_source([grid.CX(NU_ind{1}), grid.CY(NU_ind{1})], model);

  INC = zeros(2*Slength+Ulength+1, 1);

  % find maximum index:
  [maxval, maxi] = max(S);

%  disp(['maximum value: ', num2str(maxval), '  maximum cell index: ', num2str(maxi)]);
  diff_max = num_diff_flux_s.G(maxi, :);
  if any(diff_max < 0);
    disp('all values need to be positive:')
    disp(num2str(num_diff_flux_s.G(maxi, :)));
    keyboard;
  end
%  ELmax = grid.EL(maxi, :);
  Umax = U(model_data.gEI(maxi,:))'.*(grid.NX(maxi,:) + grid.NY(maxi,:));
  conv_max = num_conv_flux_s.G(maxi, :);
  Uconv_neg = zeros(size(num_conv_flux_s.G(maxi,:)));
  Uconv_pos = zeros(size(num_conv_flux_s.G(maxi,:)));
  Uconv_neg(Umax < 0) = conv_max(Umax < 0);
  fUk       = model.two_phase.waterflow([], S(maxi), model);
  Uconv_pos(Umax < 0) = - Umax(Umax < 0) .* fUk;
  if any(Uconv_pos + Uconv_neg < 0)
    disp('all values need to be positive:')
    disp(num2str(Uconv_pos + Uconv_neg));
    keyboard;
  end
  cmax = model.two_phase.injection_concentration(model);
  source_contri = (fUk - cmax) * S_upper(maxi);
  if fUk > cmax
    dodisp = source_contri < 0;
    negpos = 'positive';
  else
    dodisp = source_contri > 0;
    negpos = 'negative';
  end
  if dodisp
    disp(['source contribution should be ', negpos, ': ', num2str(source_contri)]);
    keyboard;
  end

  %sink_contri_196 = (S(196) - model.two_phase.waterflow([], S(196), model)) * S_lower(196);
%  disp(['sink at entity 196 = ', num2str(sink_contri_196)]);


  if sum(Uconv_pos + Uconv_neg) + source_contri ~= sum(conv_max) - cmax * S_upper(maxi)
%    disp('convective flux differs!');
%    keyboard;
  end

%  disp(['f(S_k) = ', num2str(fUk), ' c = ', num2str(cmax)]);

  INC(1:Slength) = grid.Ainv(NU_ind{1}, :) .* ...
    ( sum(num_diff_flux_s.G(NU_ind{1},:), 2) ...
     + sum(num_conv_flux_s.G(NU_ind{1},:), 2)  ...
    - model.two_phase.injection_concentration(model) .* S_upper ...
    + S(NU_ind{1}) .* S_lower );

  UKL = -num_diff_flux_u.G(NU_ind{2});

  INC(Slength+1:Ulength+Slength) = UKL - U;

  INC(Slength+Ulength+1:end-1) = sum(U(model_data.gEI).*(grid.NX+grid.NY), 2) ...
    - S_upper + S_lower;

%  INC(end) = sum(grid.A(NU_ind{3}, :) .* P);

  INC(end) = P(1);

  if nargout > 1

    s_diff_s = fv_frechet_operators_diff_implicit_gradient3(model, model_data, S, NU_ind{1});
    s_conv_s_u = fv_frechet_operators_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind{1}, NU_ind{2});

    u_diff_p_s = fv_frechet_operators_diff_flux_pressure_gradient(model, model_data, P, S, NU_ind{3}, NU_ind{1});


    n = 2*Slength+Ulength+1;
    m = 2*Slength+Ulength;

    sp11 = grid.Ainv(1).*sparse(s_diff_s.si, s_diff_s.sj, s_diff_s.svals, n, m);
    sp121 = sparse(s_conv_s_u.si, s_conv_s_u.sj, s_conv_s_u.svals, n, m);
    sp122 = sparse(s_conv_s_u.ui, s_conv_s_u.uj+Slength, s_conv_s_u.uvals, n, m);
    sp12 = grid.Ainv(1).*(sp121+sp122);

    sli = find(S_lower);
    sp13 = sparse(sli, sli, S_lower(sli) .* grid.Ainv(sli), n, m);

    sp1 = sp11 + sp12 + sp13;
    %sp1 = sp12 + sp13;
    sp21 = sparse(u_diff_p_s.si+Slength, u_diff_p_s.sj, u_diff_p_s.svals, n, m);
    sp22 = sparse(u_diff_p_s.pi+Slength, u_diff_p_s.pj+Slength+Ulength, u_diff_p_s.pvals, n, m);
    sp2 = -(sp21+sp22);

    sp3i = repmat(1:Slength, 1, 4)+Ulength+Slength;
    sp3j = model_data.gEI(:)+Slength;
    sp3vals = grid.NX+grid.NY;
    sp3 = sparse(sp3i,sp3j,sp3vals, n, m);

%    sp4 = sparse(2*Slength+Ulength+1*ones(1,Slength), Slength+Ulength+1:2*Slength+Ulength, grid.A(:), n, m);
    sp4 = sparse(2*Slength+Ulength+1, Slength+Ulength+1, 1, n, m);

    sp5 = sparse((1:Ulength)+Slength, (1:Ulength)+Slength, -1, n,m);

    J = sp1 + sp2 + sp3 + sp4 + sp5;
    J(isnan(J)) = 0;
  end

end