1 function diffusivity = diffusivity_buckley_leverett_simple(glob, U, params)
2 %
function diffusivity = diffusivity_buckley_leverett_simple(glob, U, params)
4 %
function computing an nonlinear diffusivity
for the buckley leverett problem
5 % ``k(u) = K p_c
'(u) \cdot \lambda_2(u) f(u) ``
7 % ``p_c(u) = u^{-\lambda},``
8 % ``p_c'(u) = -\lambda u^{-\lambda-1},``
9 % ``f(u) = \frac{\lambda_w(u)}{\lambda_w(u)+\lambda_n(u)},``
10 % ``\lambda_w(u) = \frac{u^3}{\mu_1},``
11 % ``\lambda_n(u) = \frac{(1-u)^3}{\mu_2}.``
14 % U : solution at points in
'glob'
16 % required fields of params:
17 % diff_K : diffussion factor `K`
18 % bl_mu_1 : viscosity factor `\mu_1`
19 % bl_mu_2 : viscosity factor `\mu_2`
20 % bl_lambda : mobility parameter `\lambda`
22 % See also diffusivity_buckley_leverett_simple_derivative().
24 % generated fields of diffusivity:
25 % epsilon: upper bound on diffusivity value
26 % K: vector with diffusivity values
30 if ~isempty(glob) && size(glob,1) < size(glob,2)
31 warning('coordinates in variable glob are given row-wise, but expected them to be column-wise');
37 if params.newton_regularisation && ( ( max(U) > 1 ) ...% && max(U)-1e7*eps < 1 )...
38 || ( min(U) < 0 ) )%&& min(U)+1e7*eps > 0 ) )
43 if( params.debug && any(U<0))
44 error('argument function is negative, which is disallowed here');
47 ld = params.bl_lambda;
52 lambda2 = (1-U).^3/mu2;
54 mobility = lambda1 + lambda2;
56 %deg1 = (1-U).^(-0.5);
57 %deg1(isinf(deg1)) = realmax;
59 diffusivity.K = + params.diff_K * ld/mu1 ...
61 .* (real(U.^(2-ld)) ./ mobility);
64 diffusivity.epsilon = max(diffusivity.K);
67 if ( any(isnan(diffusivity.K) | isinf(diffusivity.K)) )
68 error('diffusion factor has invalid elements');
69 elseif ( max(abs(imag(U))) > 1e-6 )
70 error('U has non-trivial imaginary addend.');
74 if params.decomp_mode>0
75 error('function is nonlinear and does not support affine decomposition!');