1 function [INC, J] = fv_two_phase_es(model, model_data, U, CS, NU_ind)
2 %
function INC = fv_two_phase_es(model, model_data, U, CS, NU_ind)
4 grid = model_data.grid;
6 if nargin <= 4 || isempty(NU_ind)
8 NU_ind{1} = (1:grid.nelements)
';
9 NU_ind{2} = (1:model_data.gn_edges)';
10 NU_ind{3} = (1:grid.nelements)
';
13 Slength = grid.nelements;
16 % model.diffusivity_ptr = @(glob, U, model) struct('K
', 1, 'epsilon
', 0);%model.two_phase.phi(glob, U, model);
17 % model.diffusivity_derivative_ptr = @(glob, U, model) struct('K
', 0, 'epsilon
', 0);%model.two_phase.phi_derivative(glob, U, model);
18 % model.laplacian_ptr = @(glob, U, model) -getfield(model.two_phase.Phi(glob, U, model), 'K
')';
19 % model.laplacian_derivative_ptr = @(glob, U, model) -getfield(model.two_phase.phi(glob, U, model),
'K');
21 model.diffusivity_ptr = @(glob, U, model) model.two_phase.phi(glob, U, model);
22 model.diffusivity_derivative_ptr = @(glob, U, model) model.two_phase.phi_derivative(glob, U, model);
23 model.laplacian_ptr = @(glob, U, descr) U;
24 model.laplacian_derivative_ptr = @(glob, U, descr) ones(length(U),1);
28 num_conv_flux_s = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind{1}, NU_ind{2});
30 % homogeneous neumann boundaries...
31 num_conv_flux_s.G(isnan(num_conv_flux_s.G)) = 0;
33 tmppar.local_mass_matrix = @fv_local_mass_matrix_rect;
35 tmppar.evaluate_basis = @fv_evaluate_basis;
37 S_upper = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.upper_source(elemin, loc, grid), 2, model_data.grid, tmppar);
38 S_lower = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.lower_source(elemin, loc, grid), 2, model_data.grid, tmppar);
42 [maxval, maxi] = max(S);
43 disp([
'maximum value: ', num2str(maxval),
' maximum cell index: ', num2str(maxi)]);
44 diff_max = num_diff_flux_s.G(maxi, :);
46 disp(
'all values need to be positive:')
47 disp(num2str(num_diff_flux_s.G(maxi, :)));
50 % ELmax = grid.EL(maxi, :);
51 Umax = U(model_data.gEI(maxi,:))'.*(grid.NX(maxi,:) + grid.NY(maxi,:));
52 conv_max = num_conv_flux_s.G(maxi, :);
53 Uconv_neg = zeros(size(num_conv_flux_s.G(maxi,:)));
54 Uconv_pos = zeros(size(num_conv_flux_s.G(maxi,:)));
55 Uconv_neg(Umax < 0) = conv_max(Umax < 0);
56 fUk = model.two_phase.waterflow([], S(maxi), model);
57 Uconv_pos(Umax < 0) = - Umax(Umax < 0) .* fUk;
58 if any(Uconv_pos + Uconv_neg < 0)
59 disp('all values need to be positive:')
60 disp(num2str(Uconv_pos + Uconv_neg));
63 cmax = model.two_phase.injection_concentration(model);
64 source_contri = (fUk - cmax) * S_upper(maxi);
66 dodisp = source_contri < 0;
69 dodisp = source_contri > 0;
73 disp(['source contribution should be ', negpos, ': ', num2str(source_contri)]);
78 if sum(Uconv_pos + Uconv_neg) + source_contri ~= sum(conv_max) - cmax * S_upper(maxi)
79 % disp('convective flux differs!');
86 INC = zeros(Slength, 1);
88 INC(1:Slength) = grid.Ainv(NU_ind{1}, :) .* ...
89 ( sum(num_diff_flux_s.G(NU_ind{1},:), 2) ...
90 + sum(num_conv_flux_s.G(NU_ind{1},:), 2) ...
91 - model.two_phase.injection_concentration(model) .* S_upper ...
92 + S(NU_ind{1}) .* S_lower );
103 sp11 = grid.Ainv(1).*sparse(s_diff_s.si, s_diff_s.sj, s_diff_s.svals, n, m);
104 sp121 = sparse(s_conv_s_u.si, s_conv_s_u.sj, s_conv_s_u.svals, n, m);
105 sp12 = grid.Ainv(1).*(sp121);
108 sp13 = sparse(sli, sli, S_lower(sli).*grid.Ainv(sli), n, m);
110 sp1 = sp11 + sp12 + sp13;
function num_flux_mat = fv_num_diff_flux_gradient(model, model_data, U, NU_ind)
computes a numerical diffusive flux for a diffusion problem
function spm = fv_frechet_operators_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind)
computes a jacobian of implicit non-linear convective contributions to time evolution matrices at a p...
function spm = fv_frechet_operators_diff_implicit_gradient3(model, model_data, U, NU_ind)
computes a jacobian of implicit non-linear diffusion contributions to time evolution matrices at a po...
function res = cubequadrature(poldeg, func, varargin)
integration of function func over reference cube == unit cube. by Gaussian quadrature exactly integra...