rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
fv_two_phase_es.m
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)
3 
4  grid = model_data.grid;
5 
6  if nargin <= 4 || isempty(NU_ind)
7  NU_ind = cell(1,3);
8  NU_ind{1} = (1:grid.nelements)';
9  NU_ind{2} = (1:model_data.gn_edges)';
10  NU_ind{3} = (1:grid.nelements)';
11  end
12 
13  Slength = grid.nelements;
14  S = CS;
15 
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');
20 
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);
25 
26  num_diff_flux_s = fv_num_diff_flux_gradient(model, model_data, S, NU_ind{1});
27 
28  num_conv_flux_s = fv_num_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind{1}, NU_ind{2});
29 
30  % homogeneous neumann boundaries...
31  num_conv_flux_s.G(isnan(num_conv_flux_s.G)) = 0;
32 
33  tmppar.local_mass_matrix = @fv_local_mass_matrix_rect;
34  tmppar.element_quadrature = @cubequadrature;
35  tmppar.evaluate_basis = @fv_evaluate_basis;
36  tmppar.pdeg = 0;
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);
39 
40 
41  %%%% debug start
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, :);
45  if any(diff_max < 0);
46  disp('all values need to be positive:')
47  disp(num2str(num_diff_flux_s.G(maxi, :)));
48  keyboard;
49  end
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));
61  keyboard;
62  end
63  cmax = model.two_phase.injection_concentration(model);
64  source_contri = (fUk - cmax) * S_upper(maxi);
65  if fUk > cmax
66  dodisp = source_contri < 0;
67  negpos = 'positive';
68  else
69  dodisp = source_contri > 0;
70  negpos = 'negative';
71  end
72  if dodisp
73  disp(['source contribution should be ', negpos, ': ', num2str(source_contri)]);
74  keyboard;
75  end
76 
77 
78  if sum(Uconv_pos + Uconv_neg) + source_contri ~= sum(conv_max) - cmax * S_upper(maxi)
79 % disp('convective flux differs!');
80 % keyboard;
81  end
82  %%%%% debug end
83 
84 
85 
86  INC = zeros(Slength, 1);
87 
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 );
93 
94  if nargout > 1
95 
96  s_diff_s = fv_frechet_operators_diff_implicit_gradient3(model, model_data, S, NU_ind{1});
97 
98  s_conv_s_u = fv_frechet_operators_conv_flux_waterflow_upwind(model, model_data, S, U, NU_ind{1}, NU_ind{2});
99 
100  n = Slength;
101  m = Slength;
102 
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);
106 
107  sli = find(S_lower);
108  sp13 = sparse(sli, sli, S_lower(sli).*grid.Ainv(sli), n, m);
109 
110  sp1 = sp11 + sp12 + sp13;
111 % sp1 = sp12 + sp13;
112 
113  J = sp1;
114  J(isnan(J)) = 0;
115  end
116 
117 end
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...