rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
diffusivity_buckley_leverett_simple.m
1 function diffusivity = diffusivity_buckley_leverett_simple(glob, U, params)
2 % function diffusivity = diffusivity_buckley_leverett_simple(glob, U, params)
3 %
4 % function computing an nonlinear diffusivity for the buckley leverett problem
5 % ``k(u) = K p_c'(u) \cdot \lambda_2(u) f(u) ``
6 % with
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}.``
12 %
13 % parameters:
14 % U : solution at points in 'glob'
15 %
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`
21 %
22 % See also diffusivity_buckley_leverett_simple_derivative().
23 %
24 % generated fields of diffusivity:
25 % epsilon: upper bound on diffusivity value
26 % K: vector with diffusivity values
27 
28 % glob column check
29 if params.debug
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');
32  if params.debug > 2
33  keyboard;
34  end
35  end
36 end
37 if params.newton_regularisation && ( ( max(U) > 1 ) ...% && max(U)-1e7*eps < 1 )...
38  || ( min(U) < 0 ) )%&& min(U)+1e7*eps > 0 ) )
39  U(U>1) = 1;
40  U(U<0) = 0;
41 end
42 
43 if( params.debug && any(U<0))
44  error('argument function is negative, which is disallowed here');
45 end
46 
47 ld = params.bl_lambda;
48 mu1 = params.bl_mu1;
49 mu2 = params.bl_mu2;
50 
51 lambda1 = U.^3/mu1;
52 lambda2 = (1-U).^3/mu2;
53 
54 mobility = lambda1 + lambda2;
55 
56 %deg1 = (1-U).^(-0.5);
57 %deg1(isinf(deg1)) = realmax;
58 
59 diffusivity.K = + params.diff_K * ld/mu1 ...
60  .* lambda2 ...
61  .* (real(U.^(2-ld)) ./ mobility);
62 
63 
64 diffusivity.epsilon = max(diffusivity.K);
65 
66 if params.debug
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.');
71  end
72 end
73 
74 if params.decomp_mode>0
75  error('function is nonlinear and does not support affine decomposition!');
76 end
77