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 "PorousFlowDarcyBase.h"
11 :
12 : #include "Assembly.h"
13 : #include "MooseMesh.h"
14 : #include "MooseVariable.h"
15 : #include "SystemBase.h"
16 :
17 : #include "libmesh/quadrature.h"
18 :
19 : InputParameters
20 9256 : PorousFlowDarcyBase::validParams()
21 : {
22 9256 : InputParameters params = Kernel::validParams();
23 18512 : params.addRequiredParam<RealVectorValue>("gravity",
24 : "Gravitational acceleration vector downwards (m/s^2)");
25 18512 : params.addRequiredParam<UserObjectName>(
26 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
27 18512 : params.addParam<unsigned>("full_upwind_threshold",
28 18512 : 5,
29 : "If, for each timestep, the number of "
30 : "upwind-downwind swaps in an element is less than "
31 : "this quantity, then full upwinding is used for that element. "
32 : "Otherwise the fallback scheme is employed.");
33 18512 : MooseEnum fallback_enum("quick harmonic", "quick");
34 18512 : params.addParam<MooseEnum>("fallback_scheme",
35 : fallback_enum,
36 : "quick: use nodal mobility without "
37 : "preserving mass. harmonic: use a "
38 : "harmonic mean of nodal mobilities "
39 : "and preserve fluid mass");
40 9256 : params.addClassDescription("Fully-upwinded advective Darcy flux");
41 9256 : return params;
42 9256 : }
43 :
44 5100 : PorousFlowDarcyBase::PorousFlowDarcyBase(const InputParameters & parameters)
45 : : Kernel(parameters),
46 5100 : _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
47 5100 : _dpermeability_dvar(
48 5100 : getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
49 10200 : _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
50 : "dPorousFlow_permeability_qp_dgradvar")),
51 5100 : _fluid_density_node(
52 5100 : getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
53 10200 : _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
54 : "dPorousFlow_fluid_phase_density_nodal_dvar")),
55 10200 : _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
56 10200 : _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
57 : "dPorousFlow_fluid_phase_density_qp_dvar")),
58 10200 : _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
59 5100 : _dfluid_viscosity_dvar(
60 5100 : getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
61 10200 : _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
62 10200 : _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
63 10200 : _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
64 : "dPorousFlow_grad_porepressure_qp_dgradvar")),
65 10200 : _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
66 : "dPorousFlow_grad_porepressure_qp_dvar")),
67 5100 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
68 5100 : _num_phases(_dictator.numPhases()),
69 10200 : _gravity(getParam<RealVectorValue>("gravity")),
70 5100 : _perm_derivs(_dictator.usePermDerivs()),
71 10200 : _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
72 10200 : _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
73 5100 : _proto_flux(_num_phases),
74 5100 : _jacobian(_num_phases),
75 5100 : _num_upwinds(),
76 5100 : _num_downwinds()
77 : {
78 : #ifdef LIBMESH_HAVE_TBB_API
79 : if (libMesh::n_threads() > 1)
80 : mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
81 : "accurately when using TBB and greater than 1 thread");
82 : #endif
83 5100 : }
84 :
85 : void
86 44791 : PorousFlowDarcyBase::timestepSetup()
87 : {
88 44791 : Kernel::timestepSetup();
89 44791 : _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
90 44791 : _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
91 44791 : }
92 :
93 : Real
94 292638888 : PorousFlowDarcyBase::darcyQp(unsigned int ph) const
95 : {
96 292638888 : return _grad_test[_i][_qp] *
97 292638888 : (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
98 : }
99 :
100 : Real
101 948858376 : PorousFlowDarcyBase::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
102 : {
103 948858376 : if (_dictator.notPorousFlowVariable(jvar))
104 : return 0.0;
105 :
106 948858376 : const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
107 :
108 : RealVectorValue deriv =
109 948858376 : _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
110 948858376 : _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
111 :
112 948858376 : deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
113 :
114 948858376 : if (_perm_derivs)
115 : {
116 6668872 : deriv += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
117 6668872 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
118 26675488 : for (const auto i : make_range(Moose::dim))
119 20006616 : deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
120 20006616 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
121 : }
122 :
123 948858376 : return _grad_test[_i][_qp] * deriv;
124 : }
125 :
126 : Real
127 0 : PorousFlowDarcyBase::computeQpResidual()
128 : {
129 0 : mooseError("PorousFlowDarcyBase: computeQpResidual called");
130 : return 0.0;
131 : }
132 :
133 : void
134 5330621 : PorousFlowDarcyBase::computeResidual()
135 : {
136 5330621 : computeResidualAndJacobian(JacRes::CALCULATE_RESIDUAL, 0.0);
137 5330621 : }
138 :
139 : void
140 0 : PorousFlowDarcyBase::computeJacobian()
141 : {
142 0 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, _var.number());
143 0 : }
144 :
145 : void
146 6434243 : PorousFlowDarcyBase::computeOffDiagJacobian(const unsigned int jvar)
147 : {
148 6434243 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, jvar);
149 6434243 : }
150 :
151 : void
152 11764864 : PorousFlowDarcyBase::computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
153 : {
154 11764864 : if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) && _dictator.notPorousFlowVariable(jvar))
155 0 : return;
156 :
157 : // The PorousFlow variable index corresponding to the variable number jvar
158 : const unsigned int pvar =
159 11764864 : ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
160 :
161 11764864 : prepareMatrixTag(_assembly, _var.number(), jvar);
162 11764864 : if ((_local_ke.n() == 0) && (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem
163 : // encountered in the
164 : // initial timestep when
165 : // use_displaced_mesh=true
166 : return;
167 :
168 : // The number of nodes in the element
169 11764864 : const unsigned int num_nodes = _test.size();
170 :
171 : // Compute the residual and jacobian without the mobility terms. Even if we are computing the
172 : // Jacobian we still need this in order to see which nodes are upwind and which are downwind.
173 26121565 : for (unsigned ph = 0; ph < _num_phases; ++ph)
174 : {
175 14356701 : _proto_flux[ph].assign(num_nodes, 0);
176 68415025 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
177 : {
178 346697212 : for (_i = 0; _i < num_nodes; ++_i)
179 292638888 : _proto_flux[ph][_i] += _JxW[_qp] * _coord[_qp] * darcyQp(ph);
180 : }
181 : }
182 :
183 : // for this element, record whether each node is "upwind" or "downwind" (or neither)
184 11764864 : const unsigned elem = _current_elem->id();
185 11764864 : if (_num_upwinds.find(elem) == _num_upwinds.end())
186 : {
187 2665294 : _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
188 2665294 : _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
189 2847538 : for (unsigned ph = 0; ph < _num_phases; ++ph)
190 : {
191 1514891 : _num_upwinds[elem][ph].assign(num_nodes, 0);
192 1514891 : _num_downwinds[elem][ph].assign(num_nodes, 0);
193 : }
194 : }
195 : // record the information once per nonlinear iteration
196 11764864 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
197 : {
198 6903932 : for (unsigned ph = 0; ph < _num_phases; ++ph)
199 : {
200 16893215 : for (unsigned nod = 0; nod < num_nodes; ++nod)
201 : {
202 13098238 : if (_proto_flux[ph][nod] > 0)
203 6337336 : _num_upwinds[elem][ph][nod]++;
204 6760902 : else if (_proto_flux[ph][nod] < 0)
205 6266114 : _num_downwinds[elem][ph][nod]++;
206 : }
207 : }
208 : }
209 :
210 : // based on _num_upwinds and _num_downwinds, calculate the maximum number
211 : // of upwind-downwind swaps that have been encountered in this timestep
212 : // for this element
213 11764864 : std::vector<unsigned> max_swaps(_num_phases, 0);
214 26121565 : for (unsigned ph = 0; ph < _num_phases; ++ph)
215 : {
216 64081777 : for (unsigned nod = 0; nod < num_nodes; ++nod)
217 49725076 : max_swaps[ph] = std::max(
218 49725076 : max_swaps[ph], std::min(_num_upwinds[elem][ph][nod], _num_downwinds[elem][ph][nod]));
219 : }
220 :
221 : // size the _jacobian correctly and calculate it for the case residual = _proto_flux
222 11764864 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
223 : {
224 14463114 : for (unsigned ph = 0; ph < _num_phases; ++ph)
225 : {
226 8028871 : _jacobian[ph].resize(_local_ke.m());
227 35726889 : for (_i = 0; _i < _test.size(); _i++)
228 : {
229 27698018 : _jacobian[ph][_i].assign(_local_ke.n(), 0.0);
230 149756918 : for (_j = 0; _j < _phi.size(); _j++)
231 1070917276 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
232 948858376 : _jacobian[ph][_i][_j] += _JxW[_qp] * _coord[_qp] * darcyQpJacobian(jvar, ph);
233 : }
234 : }
235 : }
236 :
237 : // Loop over all the phases, computing the mass flux, which
238 : // gets placed into _proto_flux, and the derivative of this
239 : // which gets placed into _jacobian
240 26121565 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
241 : {
242 14356701 : if (max_swaps[ph] < _full_upwind_threshold)
243 14355105 : fullyUpwind(res_or_jac, ph, pvar);
244 : else
245 : {
246 1596 : switch (_fallback_scheme)
247 : {
248 596 : case FallbackEnum::QUICK:
249 596 : quickUpwind(res_or_jac, ph, pvar);
250 : break;
251 1000 : case FallbackEnum::HARMONIC:
252 1000 : harmonicMean(res_or_jac, ph, pvar);
253 : break;
254 : }
255 : }
256 : }
257 :
258 : // Add results to the Residual or Jacobian
259 11764864 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
260 : {
261 5330621 : prepareVectorTag(_assembly, _var.number());
262 24550651 : for (_i = 0; _i < _test.size(); _i++)
263 41247088 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
264 22027058 : _local_re(_i) += _proto_flux[ph][_i];
265 5330621 : accumulateTaggedLocalResidual();
266 :
267 5330621 : if (_has_save_in)
268 0 : for (unsigned int i = 0; i < _save_in.size(); i++)
269 0 : _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
270 : }
271 :
272 11764864 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
273 : {
274 29851853 : for (_i = 0; _i < _test.size(); _i++)
275 132547630 : for (_j = 0; _j < _phi.size(); _j++)
276 231188920 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
277 122058900 : _local_ke(_i, _j) += _jacobian[ph][_i][_j];
278 :
279 6434243 : accumulateTaggedLocalMatrix();
280 :
281 6434243 : if (_has_diag_save_in && jvar == _var.number())
282 : {
283 : unsigned int rows = _local_ke.m();
284 0 : DenseVector<Number> diag(rows);
285 0 : for (unsigned int i = 0; i < rows; i++)
286 0 : diag(i) = _local_ke(i, i);
287 :
288 0 : for (unsigned int i = 0; i < _diag_save_in.size(); i++)
289 0 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
290 : }
291 : }
292 11764864 : }
293 :
294 : void
295 14355105 : PorousFlowDarcyBase::fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
296 : {
297 : /**
298 : * Perform the full upwinding by multiplying the residuals at the upstream nodes by their
299 : * mobilities.
300 : * Mobility is different for each phase, and in each situation:
301 : * mobility = density / viscosity for single-component Darcy flow
302 : * mobility = mass_fraction * density * relative_perm / viscosity for multi-component,
303 : *multiphase flow
304 : * mobility = enthalpy * density * relative_perm / viscosity for heat convection
305 : *
306 : * The residual for the kernel is the sum over Darcy fluxes for each phase.
307 : * The Darcy flux for a particular phase is
308 : * R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
309 : * for node i. where int is the integral over the element.
310 : * However, in fully-upwind, the first step is to take the mobility outside the integral,
311 : * which was done in the _proto_flux calculation above.
312 : *
313 : * NOTE: Physically _proto_flux[i][ph] is a measure of fluid of phase ph flowing out of node i.
314 : * If we had left in mobility, it would be exactly the component mass flux flowing out of node
315 : *i.
316 : *
317 : * This leads to the definition of upwinding:
318 : *
319 : * If _proto_flux(i)[ph] is positive then we use mobility_i. That is we use the upwind value of
320 : * mobility.
321 : *
322 : * The final subtle thing is we must also conserve fluid mass: the total component mass flowing
323 : * out of node i must be the sum of the masses flowing into the other nodes.
324 : **/
325 :
326 : // The number of nodes in the element
327 14355105 : const unsigned int num_nodes = _test.size();
328 :
329 : Real mob;
330 : Real dmob;
331 : // Define variables used to ensure mass conservation
332 : Real total_mass_out = 0.0;
333 : Real total_in = 0.0;
334 :
335 : // The following holds derivatives of these
336 : std::vector<Real> dtotal_mass_out;
337 : std::vector<Real> dtotal_in;
338 14355105 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
339 : {
340 8028643 : dtotal_mass_out.assign(num_nodes, 0.0);
341 8028643 : dtotal_in.assign(num_nodes, 0.0);
342 : }
343 :
344 : // Perform the upwinding using the mobility
345 14355105 : std::vector<bool> upwind_node(num_nodes);
346 64073989 : for (unsigned int n = 0; n < num_nodes; ++n)
347 : {
348 49718884 : if (_proto_flux[ph][n] >= 0.0) // upstream node
349 : {
350 : upwind_node[n] = true;
351 : // The mobility at the upstream node
352 25781776 : mob = mobility(n, ph);
353 25781776 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
354 : {
355 : // The derivative of the mobility wrt the PorousFlow variable
356 14410196 : dmob = dmobility(n, ph, pvar);
357 :
358 78133518 : for (_j = 0; _j < _phi.size(); _j++)
359 63723322 : _jacobian[ph][n][_j] *= mob;
360 :
361 14410196 : if (_test.size() == _phi.size())
362 : /* mobility at node=n depends only on the variables at node=n, by construction. For
363 : * linear-lagrange variables, this means that Jacobian entries involving the derivative
364 : * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
365 : * [n][n] in the line below. However, for other variable types (eg constant monomials)
366 : * I cannot tell what variable number contributes to the derivative. However, in all
367 : * cases I can possibly imagine, the derivative is zero anyway, since in the full
368 : * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
369 : */
370 12574731 : _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
371 :
372 78133518 : for (_j = 0; _j < _phi.size(); _j++)
373 63723322 : dtotal_mass_out[_j] += _jacobian[ph][n][_j];
374 : }
375 25781776 : _proto_flux[ph][n] *= mob;
376 25781776 : total_mass_out += _proto_flux[ph][n];
377 : }
378 : else
379 : {
380 : upwind_node[n] = false;
381 23937108 : total_in -= _proto_flux[ph][n]; /// note the -= means the result is positive
382 23937108 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
383 71619512 : for (_j = 0; _j < _phi.size(); _j++)
384 58332506 : dtotal_in[_j] -= _jacobian[ph][n][_j];
385 : }
386 : }
387 :
388 : // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
389 : // weighted by their proto_flux values
390 64073989 : for (unsigned int n = 0; n < num_nodes; ++n)
391 : {
392 49718884 : if (!upwind_node[n]) // downstream node
393 : {
394 23937108 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
395 71619512 : for (_j = 0; _j < _phi.size(); _j++)
396 : {
397 58332506 : _jacobian[ph][n][_j] *= total_mass_out / total_in;
398 58332506 : _jacobian[ph][n][_j] +=
399 58332506 : _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
400 58332506 : dtotal_in[_j] * total_mass_out / total_in / total_in);
401 : }
402 23937108 : _proto_flux[ph][n] *= total_mass_out / total_in;
403 : }
404 : }
405 14355105 : }
406 :
407 : void
408 596 : PorousFlowDarcyBase::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
409 : {
410 : // The number of nodes in the element
411 596 : const unsigned int num_nodes = _test.size();
412 :
413 : Real mob;
414 : Real dmob;
415 :
416 : // Use the raw nodal mobility
417 2788 : for (unsigned int n = 0; n < num_nodes; ++n)
418 : {
419 : // The mobility at the node
420 2192 : mob = mobility(n, ph);
421 2192 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
422 : {
423 : // The derivative of the mobility wrt the PorousFlow variable
424 336 : dmob = dmobility(n, ph, pvar);
425 :
426 1488 : for (_j = 0; _j < _phi.size(); _j++)
427 1152 : _jacobian[ph][n][_j] *= mob;
428 :
429 336 : if (_test.size() == _phi.size())
430 : /* mobility at node=n depends only on the variables at node=n, by construction. For
431 : * linear-lagrange variables, this means that Jacobian entries involving the derivative
432 : * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
433 : * [n][n] in the line below. However, for other variable types (eg constant monomials)
434 : * I cannot tell what variable number contributes to the derivative. However, in all
435 : * cases I can possibly imagine, the derivative is zero anyway, since in the full
436 : * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
437 : */
438 336 : _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
439 : }
440 2192 : _proto_flux[ph][n] *= mob;
441 : }
442 596 : }
443 :
444 : void
445 1000 : PorousFlowDarcyBase::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
446 : {
447 : // The number of nodes in the element
448 1000 : const unsigned int num_nodes = _test.size();
449 :
450 1000 : std::vector<Real> mob(num_nodes);
451 : unsigned num_zero = 0;
452 : unsigned zero_mobility_node = std::numeric_limits<unsigned>::max();
453 : Real harmonic_mob = 0;
454 5000 : for (unsigned n = 0; n < num_nodes; ++n)
455 : {
456 4000 : mob[n] = mobility(n, ph);
457 4000 : if (mob[n] == 0.0)
458 : {
459 : zero_mobility_node = n;
460 0 : num_zero++;
461 : }
462 : else
463 4000 : harmonic_mob += 1.0 / mob[n];
464 : }
465 1000 : if (num_zero > 0)
466 : harmonic_mob = 0.0;
467 : else
468 1000 : harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
469 :
470 : // d(harmonic_mob)/d(PorousFlow variable at node n)
471 1000 : std::vector<Real> dharmonic_mob(num_nodes, 0.0);
472 1000 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
473 : {
474 120 : const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
475 120 : if (num_zero == 0)
476 600 : for (unsigned n = 0; n < num_nodes; ++n)
477 480 : dharmonic_mob[n] = dmobility(n, ph, pvar) * harm2 / std::pow(mob[n], 2);
478 0 : else if (num_zero == 1)
479 0 : dharmonic_mob[zero_mobility_node] =
480 0 : num_nodes * dmobility(zero_mobility_node, ph, pvar); // other derivs are zero
481 : // if num_zero > 1 then all dharmonic_mob = 0.0
482 : }
483 :
484 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
485 600 : for (unsigned n = 0; n < num_nodes; ++n)
486 2400 : for (_j = 0; _j < _phi.size(); _j++)
487 : {
488 1920 : _jacobian[ph][n][_j] *= harmonic_mob;
489 1920 : if (_test.size() == _phi.size())
490 1920 : _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
491 : }
492 :
493 1000 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
494 4400 : for (unsigned n = 0; n < num_nodes; ++n)
495 3520 : _proto_flux[ph][n] *= harmonic_mob;
496 1000 : }
497 :
498 : Real
499 0 : PorousFlowDarcyBase::mobility(unsigned nodenum, unsigned phase) const
500 : {
501 0 : return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
502 : }
503 :
504 : Real
505 0 : PorousFlowDarcyBase::dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
506 : {
507 0 : Real dm = _dfluid_density_node_dvar[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
508 0 : dm -= _fluid_density_node[nodenum][phase] * _dfluid_viscosity_dvar[nodenum][phase][pvar] /
509 : std::pow(_fluid_viscosity[nodenum][phase], 2);
510 0 : return dm;
511 : }
|
