test_two_phase_datafunc.m

Go to the documentation of this file.

function OK = test_two_phase_datafunc
% performing a test of the correct implementation of derivatives in
% TwoPhaseData class.
%
% It generates a simple '2x2' on which the gradients of the functions are
% computed numerically.
%
% Return values:
%  OK: '1', if test is OK, '0' otherwise

% Martin Drohmann 21.05.2008

  OK = 1;

  tpd = TwoPhaseData.Michel;
  OK = check_two_phase(tpd, 'michel') & OK;
  tpd = TwoPhaseData.BrooksCorey;
  OK = check_two_phase(tpd, 'brooks_corey') & OK;
  tpd = TwoPhaseData.MichelSimple;
  OK = check_two_phase(tpd, 'michel_simple') & OK;
end

function OK = check_two_phase(tpd, title_prefix)
  descr = tpd.default_descr([]);
  descr.h = 5e-7;

  OK = check_convergence(@water_permeability, @water_permeability_derivative, tpd, descr, [title_prefix, '_water']);
  OK = check_convergence(@oil_permeability, @oil_permeability_derivative, tpd, descr, [title_prefix, '_oil']) & OK;
  OK = check_convergence(@capillary_pressure, @capillary_pressure_derivative, tpd, descr, [title_prefix, '_pressure']) & OK;
  OK = check_convergence(@capillary_pressure_derivative, @capillary_pressure_second_derivative, tpd, descr, [title_prefix, '_pressure_deriv']) & OK;
  OK = check_convergence(@phi, @phi_derivative, tpd, descr, [title_prefix, '_phi']) & OK;
  OK = check_convergence(@Phi, @phi, tpd, descr, [title_prefix, '_Phi']) & OK;
  OK = check_convergence(@waterflow, @waterflow_derivative, tpd, descr, [title_prefix, '_waterflow']) & OK;
  OK = check_convergence(@total_mobility, @total_mobility_deriv, tpd, descr, [title_prefix, '_mobility']) & OK;
end


function [OK] = check_convergence(f,df,tpd, descr, title)
  % EOC convergence check of discrete gradient computations
  OK = 1;

  maxerr = zeros(1,7);
  h = 0.1;

  for refstep=1:7

    c_df = f(tpd, [], 0.1:h:0.9, descr);
    if isfield(c_df, 'K')
      c_df = diff(c_df.K) / h;
    else
      c_df = diff(c_df) / h;
    end
    df_c = df(tpd, [], 0.1+h/2:h:0.9-h/2, descr);
    if isfield(df_c, 'K')
      df_c = df_c.K;
    end

    maxerr(refstep) = max(abs(df_c - c_df));

    h = h/2;

  end
  eocrefstep = log(maxerr(1:end-1)./maxerr(2:end))/log(2);
  if any(eocrefstep(1:3) ~=0) && (min(eocrefstep(3:end)) < 1.5 && max(maxerr(3:end)) > 1e-13)
    disp(['gradient for ', title, ' incorrect!']);
    disp(['EOC sequence has too small entries(<1.5): ', mat2str(eocrefstep)]);
    keyboard;
    OK = 0;
  end
end