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 4714 : PorousFlowDarcyBase::validParams()
21 : {
22 4714 : InputParameters params = Kernel::validParams();
23 9428 : params.addRequiredParam<RealVectorValue>("gravity",
24 : "Gravitational acceleration vector downwards (m/s^2)");
25 9428 : params.addRequiredParam<UserObjectName>(
26 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
27 9428 : params.addParam<unsigned>("full_upwind_threshold",
28 9428 : 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 9428 : MooseEnum fallback_enum("quick harmonic", "quick");
34 9428 : 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 4714 : params.addClassDescription("Fully-upwinded advective Darcy flux");
41 4714 : return params;
42 4714 : }
43 :
44 2537 : PorousFlowDarcyBase::PorousFlowDarcyBase(const InputParameters & parameters)
45 : : Kernel(parameters),
46 2537 : _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
47 2537 : _dpermeability_dvar(
48 2537 : getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
49 5074 : _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
50 : "dPorousFlow_permeability_qp_dgradvar")),
51 2537 : _fluid_density_node(
52 2537 : getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
53 5074 : _dfluid_density_node_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
54 : "dPorousFlow_fluid_phase_density_nodal_dvar")),
55 5074 : _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
56 5074 : _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
57 : "dPorousFlow_fluid_phase_density_qp_dvar")),
58 5074 : _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")),
59 2537 : _dfluid_viscosity_dvar(
60 2537 : getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
61 5074 : _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_nodal")),
62 5074 : _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
63 5074 : _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
64 : "dPorousFlow_grad_porepressure_qp_dgradvar")),
65 5074 : _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
66 : "dPorousFlow_grad_porepressure_qp_dvar")),
67 2537 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
68 2537 : _num_phases(_dictator.numPhases()),
69 5074 : _gravity(getParam<RealVectorValue>("gravity")),
70 2537 : _perm_derivs(_dictator.usePermDerivs()),
71 5074 : _full_upwind_threshold(getParam<unsigned>("full_upwind_threshold")),
72 5074 : _fallback_scheme(getParam<MooseEnum>("fallback_scheme").getEnum<FallbackEnum>()),
73 2537 : _proto_flux(_num_phases),
74 2537 : _jacobian(_num_phases),
75 2537 : _num_upwinds(),
76 2537 : _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 2537 : }
84 :
85 : void
86 24939 : PorousFlowDarcyBase::timestepSetup()
87 : {
88 24939 : Kernel::timestepSetup();
89 24939 : _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
90 24939 : _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
91 24939 : }
92 :
93 : Real
94 199249984 : PorousFlowDarcyBase::darcyQp(unsigned int ph) const
95 : {
96 199249984 : return _grad_test[_i][_qp] *
97 199249984 : (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
98 : }
99 :
100 : Real
101 648264144 : PorousFlowDarcyBase::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
102 : {
103 648264144 : if (_dictator.notPorousFlowVariable(jvar))
104 : return 0.0;
105 :
106 648264144 : const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
107 :
108 : RealVectorValue deriv =
109 648264144 : _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
110 648264144 : _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
111 :
112 648264144 : deriv += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
113 :
114 648264144 : if (_perm_derivs)
115 : {
116 4853504 : deriv += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
117 4853504 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
118 19414016 : for (const auto i : make_range(Moose::dim))
119 14560512 : deriv += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
120 14560512 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
121 : }
122 :
123 648264144 : 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 3542371 : PorousFlowDarcyBase::computeResidual()
135 : {
136 3542371 : computeResidualAndJacobian(JacRes::CALCULATE_RESIDUAL, 0.0);
137 3542371 : }
138 :
139 : void
140 336 : PorousFlowDarcyBase::computeJacobian()
141 : {
142 336 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, _var.number());
143 336 : }
144 :
145 : void
146 4281104 : PorousFlowDarcyBase::computeOffDiagJacobian(const unsigned int jvar)
147 : {
148 4281104 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, jvar);
149 4281104 : }
150 :
151 : void
152 7823811 : PorousFlowDarcyBase::computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
153 : {
154 7823811 : 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 7823811 : ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
160 :
161 7823811 : prepareMatrixTag(_assembly, _var.number(), jvar);
162 7823811 : 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 7823811 : 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 17402338 : for (unsigned ph = 0; ph < _num_phases; ++ph)
174 : {
175 9578527 : _proto_flux[ph].assign(num_nodes, 0);
176 46038679 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
177 : {
178 235710136 : for (_i = 0; _i < num_nodes; ++_i)
179 199249984 : _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 7823811 : const unsigned elem = _current_elem->id();
185 7823811 : if (_num_upwinds.find(elem) == _num_upwinds.end())
186 : {
187 1793238 : _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
188 1793238 : _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
189 1916089 : for (unsigned ph = 0; ph < _num_phases; ++ph)
190 : {
191 1019470 : _num_upwinds[elem][ph].assign(num_nodes, 0);
192 1019470 : _num_downwinds[elem][ph].assign(num_nodes, 0);
193 : }
194 : }
195 : // record the information once per nonlinear iteration
196 7823811 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
197 : {
198 4606328 : for (unsigned ph = 0; ph < _num_phases; ++ph)
199 : {
200 11390969 : for (unsigned nod = 0; nod < num_nodes; ++nod)
201 : {
202 8857614 : if (_proto_flux[ph][nod] > 0)
203 4291514 : _num_upwinds[elem][ph][nod]++;
204 4566100 : else if (_proto_flux[ph][nod] < 0)
205 4235652 : _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 7823811 : std::vector<unsigned> max_swaps(_num_phases, 0);
214 17402338 : for (unsigned ph = 0; ph < _num_phases; ++ph)
215 : {
216 43095427 : for (unsigned nod = 0; nod < num_nodes; ++nod)
217 33516900 : max_swaps[ph] = std::max(
218 33516900 : 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 7823811 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
223 : {
224 9633744 : for (unsigned ph = 0; ph < _num_phases; ++ph)
225 : {
226 5352304 : _jacobian[ph].resize(_local_ke.m());
227 23989652 : for (_i = 0; _i < _test.size(); _i++)
228 : {
229 18637348 : _jacobian[ph][_i].assign(_local_ke.n(), 0.0);
230 101594640 : for (_j = 0; _j < _phi.size(); _j++)
231 731221436 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
232 648264144 : _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 17402338 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
241 : {
242 9578527 : if (max_swaps[ph] < _full_upwind_threshold)
243 9577295 : fullyUpwind(res_or_jac, ph, pvar);
244 : else
245 : {
246 1232 : switch (_fallback_scheme)
247 : {
248 432 : case FallbackEnum::QUICK:
249 432 : quickUpwind(res_or_jac, ph, pvar);
250 : break;
251 800 : case FallbackEnum::HARMONIC:
252 800 : harmonicMean(res_or_jac, ph, pvar);
253 : break;
254 : }
255 : }
256 : }
257 :
258 : // Add results to the Residual or Jacobian
259 7823811 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
260 : {
261 3542371 : prepareVectorTag(_assembly, _var.number());
262 16479447 : for (_i = 0; _i < _test.size(); _i++)
263 27816628 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
264 14879552 : _local_re(_i) += _proto_flux[ph][_i];
265 3542371 : accumulateTaggedLocalResidual();
266 :
267 3542371 : 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 7823811 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
273 : {
274 20039732 : for (_i = 0; _i < _test.size(); _i++)
275 90004704 : for (_j = 0; _j < _phi.size(); _j++)
276 157203704 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
277 82957292 : _local_ke(_i, _j) += _jacobian[ph][_i][_j];
278 :
279 4281440 : accumulateTaggedLocalMatrix();
280 :
281 4281440 : 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 7823811 : }
293 :
294 : void
295 9577295 : 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 9577295 : 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 9577295 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
339 : {
340 5352144 : dtotal_mass_out.assign(num_nodes, 0.0);
341 5352144 : dtotal_in.assign(num_nodes, 0.0);
342 : }
343 :
344 : // Perform the upwinding using the mobility
345 9577295 : std::vector<bool> upwind_node(num_nodes);
346 43089331 : for (unsigned int n = 0; n < num_nodes; ++n)
347 : {
348 33512036 : if (_proto_flux[ph][n] >= 0.0) // upstream node
349 : {
350 : upwind_node[n] = true;
351 : // The mobility at the upstream node
352 17370475 : mob = mobility(n, ph);
353 17370475 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
354 : {
355 : // The derivative of the mobility wrt the PorousFlow variable
356 9698681 : dmob = dmobility(n, ph, pvar);
357 :
358 53004943 : for (_j = 0; _j < _phi.size(); _j++)
359 43306262 : _jacobian[ph][n][_j] *= mob;
360 :
361 9698681 : 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 8458137 : _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
371 :
372 53004943 : for (_j = 0; _j < _phi.size(); _j++)
373 43306262 : dtotal_mass_out[_j] += _jacobian[ph][n][_j];
374 : }
375 17370475 : _proto_flux[ph][n] *= mob;
376 17370475 : total_mass_out += _proto_flux[ph][n];
377 : }
378 : else
379 : {
380 : upwind_node[n] = false;
381 16141561 : total_in -= _proto_flux[ph][n]; /// note the -= means the result is positive
382 16141561 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
383 48586721 : for (_j = 0; _j < _phi.size(); _j++)
384 39648662 : 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 43089331 : for (unsigned int n = 0; n < num_nodes; ++n)
391 : {
392 33512036 : if (!upwind_node[n]) // downstream node
393 : {
394 16141561 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
395 48586721 : for (_j = 0; _j < _phi.size(); _j++)
396 : {
397 39648662 : _jacobian[ph][n][_j] *= total_mass_out / total_in;
398 39648662 : _jacobian[ph][n][_j] +=
399 39648662 : _proto_flux[ph][n] * (dtotal_mass_out[_j] / total_in -
400 39648662 : dtotal_in[_j] * total_mass_out / total_in / total_in);
401 : }
402 16141561 : _proto_flux[ph][n] *= total_mass_out / total_in;
403 : }
404 : }
405 9577295 : }
406 :
407 : void
408 432 : PorousFlowDarcyBase::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
409 : {
410 : // The number of nodes in the element
411 432 : const unsigned int num_nodes = _test.size();
412 :
413 : Real mob;
414 : Real dmob;
415 :
416 : // Use the raw nodal mobility
417 2096 : for (unsigned int n = 0; n < num_nodes; ++n)
418 : {
419 : // The mobility at the node
420 1664 : mob = mobility(n, ph);
421 1664 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
422 : {
423 : // The derivative of the mobility wrt the PorousFlow variable
424 224 : dmob = dmobility(n, ph, pvar);
425 :
426 1056 : for (_j = 0; _j < _phi.size(); _j++)
427 832 : _jacobian[ph][n][_j] *= mob;
428 :
429 224 : 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 224 : _jacobian[ph][n][n] += dmob * _proto_flux[ph][n];
439 : }
440 1664 : _proto_flux[ph][n] *= mob;
441 : }
442 432 : }
443 :
444 : void
445 800 : PorousFlowDarcyBase::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
446 : {
447 : // The number of nodes in the element
448 800 : const unsigned int num_nodes = _test.size();
449 :
450 800 : 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 4000 : for (unsigned n = 0; n < num_nodes; ++n)
455 : {
456 3200 : mob[n] = mobility(n, ph);
457 3200 : if (mob[n] == 0.0)
458 : {
459 : zero_mobility_node = n;
460 0 : num_zero++;
461 : }
462 : else
463 3200 : harmonic_mob += 1.0 / mob[n];
464 : }
465 800 : if (num_zero > 0)
466 : harmonic_mob = 0.0;
467 : else
468 800 : harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
469 :
470 : // d(harmonic_mob)/d(PorousFlow variable at node n)
471 800 : std::vector<Real> dharmonic_mob(num_nodes, 0.0);
472 800 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
473 : {
474 96 : const Real harm2 = std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
475 96 : if (num_zero == 0)
476 480 : for (unsigned n = 0; n < num_nodes; ++n)
477 384 : 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 480 : for (unsigned n = 0; n < num_nodes; ++n)
486 1920 : for (_j = 0; _j < _phi.size(); _j++)
487 : {
488 1536 : _jacobian[ph][n][_j] *= harmonic_mob;
489 1536 : if (_test.size() == _phi.size())
490 1536 : _jacobian[ph][n][_j] += dharmonic_mob[_j] * _proto_flux[ph][n];
491 : }
492 :
493 800 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
494 3520 : for (unsigned n = 0; n < num_nodes; ++n)
495 2816 : _proto_flux[ph][n] *= harmonic_mob;
496 800 : }
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 : }
|
