Interface.m

classdef Interface

  methods

    function diff = Phi(this, glob, S, model)

      max_S = max(S(:));
%      S_sorted = sort(S(:));

      if ~isfield(model, 'h')
        h = 1e-6; %max(min(abs(S_sorted(2:end) - S_sorted(1:end-1))), 1e-3);
      else
        h = model.h;
      end

      intrange = 0:h:max_S;
      intlength = length(intrange)-1;
      int_args = phi(this, [], intrange, model);
      if isnan(int_args.K(1))
        int_args.K(1) = 0;
      end

      ints = cumsum(int_args.K .* [1, ones(1, intlength)*h]);
      S_indices = floor(S/h)+1;

      diff.K = ints(S_indices);
      diff.epsilon = max(diff.K);
    end

    function diff = phi(this, glob, S, model)
      % function diff = phi(this, glob, S, model)
      % returns the function `\phi(s)`.
      %
      kw = water_permeability(this, glob, S, model);

      ko = oil_permeability(this, glob, S, model);

      pc_derivative = capillary_pressure_derivative(this, glob, S, model);

      diff.K = kw .* ko ./ (kw + ko) .* pc_derivative;
      diff.epsilon = max(diff.K);
    end

    function ddiff = phi_derivative(this, glob, S, model)
      % function ddiff = phi_derivative(this, glob, S, model)
      % returns the derivative `\phi'(s)`.
      kw = water_permeability(this, glob, S, model);

      dkw = water_permeability_derivative(this, glob, S, model);

      ko = oil_permeability(this, glob, S, model);

      dko = oil_permeability_derivative(this, glob, S, model);

      pc_derivative = capillary_pressure_derivative(this, glob, S, model);

      pc_sec_derivative = capillary_pressure_second_derivative(this, glob, S, model);

      oben = ko .* kw .* pc_derivative;

      oben_der = dko .* kw .* pc_derivative + ko .* dkw .* pc_derivative + ko .* kw .* pc_sec_derivative;

      unten = kw + ko;

      unten_der = dkw + dko;

      nominator = (oben_der .* unten - unten_der .* oben);
      denominator = unten .* unten;
      denominator(nominator == 0) = 1;
      ddiff.K = nominator ./ denominator;
      assert(all(~isinf(ddiff.K)));
      assert(all(~isnan(ddiff.K)));
      ddiff.epsilon = max(ddiff.K);
    end

    function fw = waterflow(this, glob, S, model)
      % function fw = waterflow(this, glob, S, model)

      kw = water_permeability(this, glob, S, model);

      ko = oil_permeability(this, glob, S, model);

      fw = kw./(kw + ko);
    end

    function dfw = waterflow_derivative(this, glob, S, model)
      % function dfw = waterflow_derivative(this, glob, S, model)

      kw        = water_permeability(this, glob, S, model);

      dkw       = water_permeability_derivative(this, glob, S, model);

      ko        = oil_permeability(this, glob, S, model);

      dko       = oil_permeability_derivative(this, glob, S, model);

      unten     = kw + ko;

      unten_der = dkw + dko;

      dfw = (dkw .* unten - unten_der .* kw) ./ (unten .* unten);

    end

    function fwU = waterflow_flux(this, glob, S, U, model)
      % function fwU = waterflow_flux(this, glob, S, U, model)
      % returns the waterflow flux over the boundaries `f_w(s) \cdot \underscore{u}`.

      kw = water_permeability(this, glob, S, model);

      assert(size(U) == [size(kw, 1), 2]);
      ko = oil_permeability(this, glob, S, model);
      scalar = kw./(kw + ko);

      fwU = [ scalar, scalar ] .* U;

    end

    function fw_deriv_U = waterflow_deriv_flux(this, glob, S, U, model)
      % function fwU = waterflow_deriv_flux(this, glob, S, model)
      % returns the waterflow flux over the boundaries, actually only `\underscore{u}`.

      fw_deriv_U = U;

    end

    function m = total_mobility(this, glob, S, model)
      % function m = total_mobility(this, glob, S, model)

      kw = water_permeability(this, glob, S, model);
      ko = oil_permeability(this, glob, S, model);
      m.K = kw + ko;
      m.epsilon = max(m.K);

    end

    function m = total_mobility_deriv(this, glob, S, model)
      % function m = total_mobility(this, glob, S, model)

      dkw = water_permeability_derivative(this, glob, S, model);
      dko = oil_permeability_derivative(this, glob, S, model);
      m.K = dkw + dko;
      m.epsilon = max(m.K);

    end


    function m = total_mobility_dual(this, glob, S, NBI, Sdir, dir_NBI, model)
      % function m = total_mobility_dual(this, glob, S, NBI, Sdir, dir_NBI, model)

      mob  = total_mobility(this, glob, S, model);
      mobn = mob(NBI);
      mob  = repmat(mob, 4, 1);

      assert(size(mob) == size(mobn));

      m = 1 ./ ( (2 ./ mob) + (2 ./ mobn) );
    end

    function p = visualize(this)
      Stest = 0:0.01:1;

      descr = [];
      descr = default_descr(this, descr);

      glob = [];

      kw  = water_permeability(this, glob, Stest, descr);
      ko  = oil_permeability(this, glob, Stest, descr);
      pc  = capillary_pressure(this, glob, Stest, descr);
      dpc = capillary_pressure_derivative(this, glob, Stest, descr);
      ph  = phi(this, glob, Stest, descr);
      Ph  = Phi(this, glob, Stest, descr);

      fw  = waterflow(this, glob, Stest, descr);
      m   = total_mobility(this, glob, Stest, descr);

      p = figure;

      subplot(2,3,1);    plot(Stest, kw, Stest, ko);     title('kw, ko');
      sp=subplot(2,3,2); plot(Stest, pc);                set(sp, 'YScale', 'log'); title('p_c');
      sp=subplot(2,3,3); plot(Stest, dpc);               set(sp, 'YScale', 'log'); title('p_c''');
      subplot(2,3,4);    plot(Stest, ph.K, Stest, Ph.K); title('\phi, \Phi');
      subplot(2,3,5);    plot(Stest, fw);                title('f');
      subplot(2,3,6);    plot(Stest, m.K);               title('M');
    end

  end

  methods (Abstract)

    kw = water_permeability(this, glob, S, model);

    dkw = water_permeability_derivative(this, glob, S, model);

    ko = oil_permeability(this, glob, S, model);

    dko = oil_permeability_derivative(this, glob, S, model);

    pc = capillary_pressure(this, glob, S, model);

    dpc = capillary_pressure_derivative(this, glob, S, model);

    ddpc = capillary_pressure_second_derivative(this, glob, S, model);

    c = injection_concentration(this, glob, model);

    s_under = lower_source(this, elemin, loc, grid);

    s_above = upper_source(this, elemin, loc, grid);

    descr = default_descr(this, descr);

  end


end