diffusivity_buckley_leverett.m

function diffusivity = diffusivity_buckley_leverett(glob, U, params)
% function diffusivity = diffusivity_buckley_leverett(glob, U, params)
%
% function computing an nonlinear diffusivity for the buckley leverett problem
% ``k(u) = K p_c'(u) \cdot \lambda_2(u) f(u) ``
% with
% ``p_c(u) = u^{-\lambda},``
% ``p_c'(u) = -\lambda u^{-\lambda-1},``
% ``f(u) = \frac{\lambda_w(u)}{\lambda_w(u)+\lambda_n(u)},``
% ``\lambda_w(u) = \frac{u^{\frac{2+3\lambda}{\lambda}}}{\mu_1},``
% ``\lambda_n(u) = \frac{(1-u)^2 (1-u^{\frac{2+\lambda}{\lambda}})}{\mu_2}.``
%
% parameters:
%      U                      : solution at points in 'glob'
%
% required fields of params:
%      diff_K                 : diffussion factor `K`
%      bl_mu_1                : viscosity factor `\mu_1`
%      bl_mu_2                : viscosity factor `\mu_2`
%      bl_lambda              : mobility parameter `\lambda`
%
% See also diffusivity_buckley_leverett_derivative().
%
% generated fields of diffusivity:
%     epsilon: upper bound on diffusivity value
%     K: vector with diffusivity values

% glob column check
if params.debug
  if ~isempty(glob) && size(glob,1) < size(glob,2)
    warning('coordinates in variable glob are given row-wise, but expected them to be column-wise');
    if params.debug > 2
      keyboard;
    end
  end
  if(any(U<0))
    error('argument function is negative, which is disallowed here');
  end
end

ld  = params.bl_lambda;
mu1 = params.bl_mu1;
mu2 = params.bl_mu2;

lambda1 = U.^((2+3*ld)/ld)/mu1;
lambda2 = (1-U).^2.*(1-U.^(2/ld+1))/mu2;

mobility = lambda1 + lambda2;

%deg1 = (1-U).^(-0.5);
%deg1(isinf(deg1)) = realmax;

diffusivity.K = + params.diff_K * ld/mu1 ...
                  .* lambda2 ...
                  .* (U.^(2-ld+2/ld) ./ mobility);

diffusivity.epsilon = max(diffusivity.K);

if params.debug
  if ( any(isnan(diffusivity.K) | isinf(diffusivity.K)) )
    error('diffusion factor has invalid elements');
  elseif ( max(abs(imag(U))) > 1e-6 )
    error('U has non-trivial imaginary addend.');
  end
end

if params.decomp_mode>0
  error('function is nonlinear and does not support affine decomposition!');
end