test_twophase_jacobian.m

Go to the documentation of this file.

function OK = test_twophase_jacobian
% 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.two_phase = tpd;
  descr.mean_ptr = @(X) harmmean(max(X, 0), 2);
  descr.mean_deriv1_ptr = @(X) 2./( (1./X(:,1) + 1./(X(:,2))).^2 .* X(:,1).^2 );
  descr.mean_deriv2_ptr = @(X) 2./( (1./X(:,1) + 1./(X(:,2))).^2 .* X(:,2).^2 );

  descr.gridtype = 'rectgrid';
  descr.xnumintervals = 3;
  descr.ynumintervals = 3;
  descr.xrange = [0,1];
  descr.yrange = [0,1];
  descr.bnd_rect_corner1 = [-1; -1];
  descr.bnd_rect_corner2 = [ 2; 2];
  descr.bnd_rect_index   = (-2);

  descr.rb_problem_type = 'TwoPhaseFlow';
  descr.debug = 0;
  descr.verbose = 1;
  descr.mu_names = {'mu'};
  descr.mu_ranges = [0 1];

  descr.newton_steps = 10;
  descr.newton_epsilon = 1e-8;
  descr.t = 1;
  descr.tstep = 1;

  dmodel = gen_detailed_model(descr);
  descr = dmodel.descr;
  model_data = gen_model_data(dmodel);

  grid = model_data.grid;
  Slength = grid.nelements;
  Ulength = model_data.gn_edges;
  Stest = (0.3:0.4/(Slength-1):0.7)';
  Ptest = (-1:2/(Slength-1):1)';
  Utest = (0.3:0.4/(Ulength-1):0.7)';
  Utest(1:model_data.gn_inner_edges) = 0.5;

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

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

  % sp11
  s_diff_s = fv_frechet_operators_diff_implicit_gradient3(descr, model_data, Stest);
  sp11 = sparse(s_diff_s.si, s_diff_s.sj, s_diff_s.svals, Slength, Slength);
  diff_flux_fun = @(model, model_data, S) sum(getfield(fv_num_diff_flux_gradient(model, model_data, S), 'G'), 2);
  OK = check_jacobian(descr, model_data, Stest, diff_flux_fun, sp11, [title_prefix, '_s_diff_s']);

  s_conv_s_u_func = @(model, model_data, S) s_conv_s_u_wrapper(model, model_data, S, Utest);

  s_conv_s_u = fv_frechet_operators_conv_flux_waterflow_upwind(descr, model_data, Stest, Utest);
  % sp121
  sp121 = sparse(s_conv_s_u.si, s_conv_s_u.sj, s_conv_s_u.svals, Slength, Slength);
  OK = check_jacobian(descr, model_data, Stest, s_conv_s_u_func, sp121, [title_prefix, '_s_conv_s_u']) & OK;

  % sp122
  u_conv_s_u_func = @(model, model_data, U) s_conv_s_u_wrapper(model, model_data, Stest, U);
  sp122 = sparse(s_conv_s_u.ui, s_conv_s_u.uj, s_conv_s_u.uvals, Slength, Ulength);
  OK = check_jacobian(descr, model_data, Utest, u_conv_s_u_func, sp122, [title_prefix, '_u_conv_s_u']) & OK;

  % sp1
  sat_L = Fv.TwoPhase.SaturationSpace;
  sp1_fun = @(model, model_data, CU) sat_L.apply(model, model_data, ...
    struct('S', CU(1:length(Stest)), 'U', CU(length(Stest)+1:end)));
  [~,~,sp1] = sp1_fun(descr, model_data, [Stest;Utest]);
  OK = check_jacobian(descr, model_data, [Stest;Utest], sp1_fun, sp1, [title_prefix, '_sp1']) & OK;

  % sp21
  u_diff_p_s = fv_frechet_operators_diff_flux_pressure_gradient(descr, model_data, Ptest, Stest);
  s_conv_p_s_func = @(model, model_data, S) getfield(fv_num_diff_flux_pressure_gradient(model, model_data, Ptest, S), 'G');
  sp21 = sparse(u_diff_p_s.si, u_diff_p_s.sj, u_diff_p_s.svals, Ulength, Slength);
  OK = check_jacobian(descr, model_data, Stest, s_conv_p_s_func, sp21, [title_prefix, '_s_conv_p_s']) & OK;

  % sp22
  p_conv_p_s_func = @(model, model_data, P) getfield(fv_num_diff_flux_pressure_gradient(model, model_data, P, Stest), 'G');
  sp22 = sparse(u_diff_p_s.pi, u_diff_p_s.pj, u_diff_p_s.pvals, Ulength, Slength);
  OK = check_jacobian(descr, model_data, Ptest, p_conv_p_s_func, sp22, [title_prefix, '_p_conv_p_s']) & OK;

  % sp2
  sat_V = Fv.TwoPhase.VelocitySpace;
  sp2_fun = @(model, model_data, CU) sat_V.apply(model, model_data, ...
    struct('S', CU(1:length(Stest)), 'P', CU(length(Stest)+1:end)));
  [~,~,sp2] = sp2_fun(descr, model_data, [Stest;Ptest]);
  OK = check_jacobian(descr, model_data, [Stest;Ptest], sp2_fun, sp2, [title_prefix, '_sp2']) & OK;

  % sp3
  sat_V = Fv.TwoPhase.DivergenceSpace;
  sp3_fun = @(model, model_data, CU) sat_V.apply(model, model_data, struct('U', CU));
  [~,~,sp3] = sp3_fun(descr, model_data, Utest);
  OK = check_jacobian(descr, model_data, Utest, sp3_fun, sp3, [title_prefix, '_sp3']) & OK;

  % sp4
  sp4_func = @(model, model_data, P) sum(grid.A .* P);
  sp4 = sparse(ones(1,Slength), 1:Slength, grid.A(:)');
  OK = check_jacobian(descr, model_data, Ptest, sp4_func, sp4, [title_prefix, '_sp4']) & OK;

end


function [OK, ok_mat] = check_jacobian(model, model_data, Utest, fun, jac_test, title, epsilon)

  if nargin <= 6
    epsilon = 5e-2;
  end
  jac_comp = zeros(size(jac_test));

  for j = 1:size(jac_test, 2);
    addend = zeros(size(Utest));
    addend(j) = 1e-6;
    jac_comp(:,j) = 1/(2e-6) .* (fun(model, model_data, Utest + addend) - fun(model, model_data, Utest - addend));
  end

  nz_elements = abs(jac_comp) > 1e6*eps;
  ok_mat = zeros(size(jac_test));
  ok_mat(nz_elements) = abs(jac_comp(nz_elements) - jac_test(nz_elements))./(abs(jac_comp(nz_elements))) > epsilon;
  ok_mat = ok_mat + (nz_elements) - (abs(full(jac_test)) > 100*eps);
  OK = all(~ok_mat(:));

  if ~OK
    disp(['Jacobian incorrect: ', title]);
    keyboard;
  end
end

function res = s_conv_s_u_wrapper(model, model_data, S, U)
  flux = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U);
  flux.G(isnan(flux.G)) = 0;
  res = sum(flux.G, 2);
end