Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Q2PSaturationFlux.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "MooseVariable.h"
15 : #include "SystemBase.h"
16 :
17 : #include "libmesh/quadrature.h"
18 :
19 : // C++ includes
20 : #include <iostream>
21 :
22 : registerMooseObject("RichardsApp", Q2PSaturationFlux);
23 :
24 : InputParameters
25 28 : Q2PSaturationFlux::validParams()
26 : {
27 28 : InputParameters params = Kernel::validParams();
28 56 : params.addRequiredParam<UserObjectName>(
29 : "fluid_density",
30 : "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
31 56 : params.addRequiredCoupledVar("porepressure_variable",
32 : "The variable representing the porepressure");
33 56 : params.addRequiredParam<UserObjectName>(
34 : "fluid_relperm",
35 : "A RichardsRelPerm UserObject (eg RichardsRelPermPowerGas) that defines the "
36 : "fluid relative permeability as a function of the saturation Variable.");
37 56 : params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
38 28 : params.addClassDescription(
39 : "Flux according to Darcy-Richards flow. The Variable of this Kernel must be the saturation");
40 28 : return params;
41 0 : }
42 :
43 14 : Q2PSaturationFlux::Q2PSaturationFlux(const InputParameters & parameters)
44 : : Kernel(parameters),
45 14 : _density(getUserObject<RichardsDensity>("fluid_density")),
46 14 : _pp(coupledValue("porepressure_variable")),
47 14 : _grad_pp(coupledGradient("porepressure_variable")),
48 14 : _pp_nodal(coupledDofValues("porepressure_variable")),
49 14 : _pp_var(coupled("porepressure_variable")),
50 14 : _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
51 28 : _viscosity(getParam<Real>("fluid_viscosity")),
52 28 : _gravity(getMaterialProperty<RealVectorValue>("gravity")),
53 28 : _permeability(getMaterialProperty<RealTensorValue>("permeability")),
54 14 : _num_nodes(0),
55 14 : _mobility(0),
56 14 : _dmobility_dp(0),
57 28 : _dmobility_ds(0)
58 : {
59 14 : }
60 :
61 : void
62 12006 : Q2PSaturationFlux::prepareNodalValues()
63 : {
64 12006 : _num_nodes = _pp_nodal.size();
65 :
66 : Real density;
67 : Real ddensity_dp;
68 : Real relperm;
69 : Real drelperm_ds;
70 :
71 12006 : _mobility.resize(_num_nodes);
72 12006 : _dmobility_dp.resize(_num_nodes);
73 12006 : _dmobility_ds.resize(_num_nodes);
74 :
75 40398 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
76 : {
77 28392 : density = _density.density(_pp_nodal[nodenum]); // fluid density at the node
78 28392 : ddensity_dp = _density.ddensity(_pp_nodal[nodenum]); // d(fluid density)/dP at the node
79 28392 : relperm = _relperm.relperm(
80 28392 : _var.dofValues()[nodenum]); // relative permeability of the fluid at node nodenum
81 28392 : drelperm_ds = _relperm.drelperm(_var.dofValues()[nodenum]); // d(relperm)/dsat
82 :
83 : // calculate the mobility and its derivatives wrt P and S
84 28392 : _mobility[nodenum] = density * relperm / _viscosity;
85 28392 : _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
86 28392 : _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
87 : }
88 12006 : }
89 :
90 : Real
91 88032 : Q2PSaturationFlux::computeQpResidual()
92 : {
93 : // note this is not the complete residual:
94 : // the upwind mobility parts get added in computeResidual
95 88032 : return _grad_test[_i][_qp] *
96 88032 : (_permeability[_qp] * (_grad_pp[_qp] - _density.density(_pp[_qp]) * _gravity[_qp]));
97 : }
98 :
99 : void
100 5535 : Q2PSaturationFlux::computeResidual()
101 : {
102 5535 : upwind(true, false, 0);
103 5535 : }
104 :
105 : void
106 33 : Q2PSaturationFlux::computeJacobian()
107 : {
108 33 : upwind(false, true, _var.number());
109 33 : }
110 :
111 : void
112 6439 : Q2PSaturationFlux::computeOffDiagJacobian(const unsigned int jvar)
113 : {
114 6439 : upwind(false, true, jvar);
115 6439 : }
116 :
117 : Real
118 211648 : Q2PSaturationFlux::computeQpJac(unsigned int dvar)
119 : {
120 : // this is just the derivative of the flux WITHOUT the upstream mobility terms
121 : // Those terms get added in during computeJacobian()
122 211648 : if (dvar == _pp_var)
123 104768 : return _grad_test[_i][_qp] *
124 104768 : (_permeability[_qp] *
125 104768 : (_grad_phi[_j][_qp] - _density.ddensity(_pp[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
126 : else
127 : return 0;
128 : }
129 :
130 : void
131 12007 : Q2PSaturationFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
132 : {
133 12007 : if (compute_jac && !(jvar == _var.number() || jvar == _pp_var))
134 1 : return;
135 :
136 : // calculate the mobility values and their derivatives
137 12006 : prepareNodalValues();
138 :
139 : // compute the residual without the mobility terms
140 : // Even if we are computing the jacobian we still need this
141 : // in order to see which nodes are upwind and which are downwind
142 12006 : prepareVectorTag(_assembly, _var.number());
143 :
144 40398 : for (_i = 0; _i < _test.size(); _i++)
145 116424 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
146 88032 : _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
147 :
148 12006 : if (compute_jac)
149 : {
150 6471 : prepareMatrixTag(_assembly, _var.number(), jvar);
151 21547 : for (_i = 0; _i < _test.size(); _i++)
152 60916 : for (_j = 0; _j < _phi.size(); _j++)
153 257488 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
154 211648 : _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
155 : }
156 :
157 : // Now perform the upwinding by multiplying the residuals at the
158 : // upstream nodes by their mobilities
159 : //
160 : // The residual for the kernel is the darcy flux.
161 : // This is
162 : // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
163 : // for node i. where int is the integral over the element.
164 : // However, in fully-upwind, the first step is to take the mobility outside the
165 : // integral, which was done in the _local_re calculation above.
166 : //
167 : // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
168 : // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
169 : //
170 : // This leads to the definition of upwinding:
171 : // ***
172 : // If _local_re(i) is positive then we use mobility_i. That is
173 : // we use the upwind value of mobility.
174 : // ***
175 : //
176 : // The final subtle thing is we must also conserve fluid mass: the total mass
177 : // flowing out of node i must be the sum of the masses flowing
178 : // into the other nodes.
179 :
180 : // FIRST:
181 : // this is a dirty way of getting around precision loss problems
182 : // and problems at steadystate where upwinding oscillates from
183 : // node-to-node causing nonconvergence.
184 : // I'm not sure if i actually need to do this in moose. Certainly
185 : // in cosflow it is necessary.
186 : // I will code a better algorithm if necessary
187 : bool reached_steady = true;
188 27100 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
189 : {
190 26608 : if (_local_re(nodenum) >= 1E-20)
191 : {
192 : reached_steady = false;
193 : break;
194 : }
195 : }
196 : reached_steady = false;
197 :
198 : // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
199 : // total mass out - used for mass conservation
200 : Real total_mass_out = 0;
201 : // total flux in
202 : Real total_in = 0;
203 :
204 : // the following holds derivatives of these
205 : std::vector<Real> dtotal_mass_out;
206 : std::vector<Real> dtotal_in;
207 12006 : if (compute_jac)
208 : {
209 6471 : dtotal_mass_out.assign(_num_nodes, 0);
210 6471 : dtotal_in.assign(_num_nodes, 0);
211 : }
212 :
213 : // PERFORM THE UPWINDING!
214 40398 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
215 : {
216 28392 : if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
217 : {
218 14326 : if (compute_jac)
219 : {
220 31062 : for (_j = 0; _j < _phi.size(); _j++)
221 23448 : _local_ke(nodenum, _j) *= _mobility[nodenum];
222 7614 : if (jvar == _var.number())
223 : // deriv wrt S
224 3842 : _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
225 : else
226 : // deriv wrt P
227 3772 : _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
228 31062 : for (_j = 0; _j < _phi.size(); _j++)
229 23448 : dtotal_mass_out[_j] += _local_ke(nodenum, _j);
230 : }
231 14326 : _local_re(nodenum) *= _mobility[nodenum];
232 14326 : total_mass_out += _local_re(nodenum);
233 : }
234 : else
235 : {
236 14066 : total_in -= _local_re(nodenum); // note the -= means the result is positive
237 14066 : if (compute_jac)
238 29854 : for (_j = 0; _j < _phi.size(); _j++)
239 22392 : dtotal_in[_j] -= _local_ke(nodenum, _j);
240 : }
241 : }
242 :
243 : // CONSERVE MASS
244 : // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
245 : if (!reached_steady)
246 40398 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
247 28392 : if (_local_re(nodenum) < 0)
248 : {
249 14066 : if (compute_jac)
250 29854 : for (_j = 0; _j < _phi.size(); _j++)
251 : {
252 22392 : _local_ke(nodenum, _j) *= total_mass_out / total_in;
253 22392 : _local_ke(nodenum, _j) +=
254 22392 : _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
255 22392 : dtotal_in[_j] * total_mass_out / total_in / total_in);
256 : }
257 14066 : _local_re(nodenum) *= total_mass_out / total_in;
258 : }
259 :
260 : // ADD RESULTS TO RESIDUAL OR JACOBIAN
261 12006 : if (compute_res)
262 : {
263 5535 : accumulateTaggedLocalResidual();
264 :
265 5535 : if (_has_save_in)
266 602 : for (unsigned int i = 0; i < _save_in.size(); i++)
267 301 : _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
268 : }
269 :
270 12006 : if (compute_jac)
271 : {
272 6471 : accumulateTaggedLocalMatrix();
273 6471 : if (_has_diag_save_in && jvar == _var.number())
274 : {
275 : const unsigned int rows = _local_ke.m();
276 103 : DenseVector<Number> diag(rows);
277 515 : for (unsigned int i = 0; i < rows; i++)
278 412 : diag(i) = _local_ke(i, i);
279 :
280 206 : for (unsigned int i = 0; i < _diag_save_in.size(); i++)
281 103 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
282 : }
283 : }
284 12006 : }
|
