SaturationSpace.m

classdef SaturationSpace < ILocalizedOperator

  properties (Constant)
    id = 'L_saturation';

    ret_vars = {'saturation'};
    arg_vars = {'saturation', 'velocity'};
    arg_vars_short = {'S', 'U'};
  end

  methods
    function [INC, INCoff, J] = apply(this, model, model_data, sim_data, NU_ind)
    %function INC = apply(this, model, model_data, sim_data, NU_ind)

      grid = model_data.grid;

      S = sim_data.S;
      U = sim_data.U;

      if nargin <= 4 || isempty(NU_ind)
        NU_ind = (1:length(S));
      end

      if ~isfield(sim_data, 'm')
        sim_data.m = length(S) + length(U);
      end

      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);

      num_conv_flux_s = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U);

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

      if model.debug > 0
        % 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)]);
      end

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

      INCoff = [];

      if nargout > 2

        s_diff_s = fv_frechet_operators_diff_implicit_gradient3(model, model_data, S, NU_ind);
        s_conv_s_u = fv_frechet_operators_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind);

        n = length(NU_ind);
        UN_ind = zeros(length(S),1);
        UN_ind(NU_ind) = 1:length(NU_ind);
        if ~isfield(sim_data, 'Uoff')
          sim_data.Uoff = length(S);
        end
        m = sim_data.m;

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

        sli = intersect(find(model_data.opdata.S_lower), NU_ind);
        sp13 = sparse(UN_ind(sli), sli, model_data.opdata.S_lower(sli) .* grid.Ainv(sli), n, m);

        J = sp11 + sp12 + sp13;
      end

    end

    function ret_size = ret_size(this, model_data, NU_ind)
      if nargin == 2 || isempty(NU_ind)
        ret_size = model_data.grid.nelements;
      else
        ret_size = length(NU_ind);
      end
    end

    function arg_size = arg_size(this, model_data)
      arg_size = 2 * model_data.grid.nelements;
    end

    function ipm = inner_product_matrix(this, model_data)
      ipm = model_data.W;
    end

    function [nmd, eind, eind_local] = compute_TM_global(this, model_data, TM)
      eind = ILocalizedOperator.compute_TM_global_edge(model_data, TM);
      grid = model_data.grid;
      glob2loc       = zeros(max(TM),1);
      glob2loc(eind) = 1:length(eind);
      eind_local     = glob2loc(TM);

      if isa(grid,'rectgrid') || isa(grid,'triagrid') || isa(grid,'onedgrid')
        nmd.grid_local_ext = gridpart(grid,eind);

        gle = nmd.grid_local_ext;
        mask = zeros(size(gle.NBI));
        mask(eind_local, :) = 1;
        [elnb_i, elnb_j] = find(gle.NBI(eind_local, :) > 0);
        elnb = sub2ind(size(mask), eind_local(elnb_i), elnb_j);
        mask(sub2ind(size(mask), gle.NBI(elnb), gle.INB(elnb))) = 1;
        gle.NBI(~mask) = -10;

%        eind_local_inv = setdiff(1:nmd.grid_local_ext.nelements, eind_local);
%        [tmp_i, tmp_j] = find(nmd.grid_local_ext.NBI(eind_local_inv, :) < 0);
%        nmd.grid_local_ext.NBI(sub2ind(size(nmd.grid_local_ext.NBI), eind_local_inv(tmp_i), tmp_j')) = -10;

        [nmd.gEI, nmd.gn_inner_edges, nmd.gn_boundary_edges] ...
          = compute_edge_indices(nmd.grid_local_ext);
        nmd.gn_edges = nmd.gn_inner_edges + nmd.gn_boundary_edges;
      else
        error('unsupported grid');
      end

      nmd.TM_global_args{1} = eind;
      vel_global = zeros(1, model_data.gn_edges);
      vel_global(model_data.gEI(TM, :)) = 1;
      nmd.TM_global_args{2} = find(vel_global);
    end

    function opdata = fill_opdata(this, rmodel, detailed_data)
      if isa(rmodel.M, 'DataTree.IdMapNode')
        MM = get_by_description(rmodel.M, 'L_saturation');
      else
        MM = rmodel.M;
      end
      opdata.S_lower = detailed_data.model_data.opdata.S_lower(detailed_data.TM(1:MM));
      opdata.S_upper = detailed_data.model_data.opdata.S_upper(detailed_data.TM(1:MM));
    end

    function opdata = copy_extract_opdata(this, opdata_copy, rmodel, Mid)
      MM = get_by_description(rmodel.M, Mid);
      opdata.S_lower = opdata_copy.S_lower(1:MM);
      opdata.S_upper = opdata_copy.S_upper(1:MM);
    end


  end
end