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 : template <bool is_ad>
20 : InputParameters
21 4833 : PorousFlowDarcyBaseTempl<is_ad>::validParams()
22 : {
23 : InputParameters params = GenericKernel<is_ad>::validParams();
24 9666 : params.addRequiredParam<RealVectorValue>("gravity",
25 : "Gravitational acceleration vector downwards (m/s^2)");
26 9666 : params.addRequiredParam<UserObjectName>(
27 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
28 9666 : params.addParam<unsigned>("full_upwind_threshold",
29 9666 : 5,
30 : "If, for each timestep, the number of "
31 : "upwind-downwind swaps in an element is less than "
32 : "this quantity, then full upwinding is used for that element. "
33 : "Otherwise the fallback scheme is employed.");
34 9666 : MooseEnum fallback_enum("quick harmonic", "quick");
35 9666 : params.addParam<MooseEnum>("fallback_scheme",
36 : fallback_enum,
37 : "quick: use nodal mobility without "
38 : "preserving mass. harmonic: use a "
39 : "harmonic mean of nodal mobilities "
40 : "and preserve fluid mass");
41 4833 : params.addClassDescription("Fully-upwinded advective Darcy flux");
42 4833 : return params;
43 4833 : }
44 :
45 : template <bool is_ad>
46 2602 : PorousFlowDarcyBaseTempl<is_ad>::PorousFlowDarcyBaseTempl(const InputParameters & parameters)
47 : : GenericKernel<is_ad>(parameters),
48 2602 : _permeability(this->template getGenericMaterialProperty<RealTensorValue, is_ad>(
49 : "PorousFlow_permeability_qp")),
50 2602 : _dpermeability_dvar(is_ad ? nullptr
51 2539 : : &this->template getMaterialProperty<std::vector<RealTensorValue>>(
52 : "dPorousFlow_permeability_qp_dvar")),
53 2602 : _dpermeability_dgradvar(
54 : is_ad ? nullptr
55 2539 : : &this->template getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
56 : "dPorousFlow_permeability_qp_dgradvar")),
57 5204 : _fluid_density_node(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
58 : "PorousFlow_fluid_phase_density_nodal")),
59 2602 : _dfluid_density_node_dvar(
60 : is_ad ? nullptr
61 2539 : : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
62 : "dPorousFlow_fluid_phase_density_nodal_dvar")),
63 5204 : _fluid_density_qp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
64 : "PorousFlow_fluid_phase_density_qp")),
65 2602 : _dfluid_density_qp_dvar(
66 : is_ad ? nullptr
67 2539 : : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
68 : "dPorousFlow_fluid_phase_density_qp_dvar")),
69 5204 : _fluid_viscosity(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
70 : "PorousFlow_viscosity_nodal")),
71 2602 : _dfluid_viscosity_dvar(
72 : is_ad ? nullptr
73 2539 : : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
74 : "dPorousFlow_viscosity_nodal_dvar")),
75 5204 : _pp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
76 : "PorousFlow_porepressure_nodal")),
77 5204 : _grad_p(this->template getGenericMaterialProperty<std::vector<RealGradient>, is_ad>(
78 : "PorousFlow_grad_porepressure_qp")),
79 2602 : _dgrad_p_dgrad_var(is_ad ? nullptr
80 2539 : : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
81 : "dPorousFlow_grad_porepressure_qp_dgradvar")),
82 2602 : _dgrad_p_dvar(is_ad
83 : ? nullptr
84 2539 : : &this->template getMaterialProperty<std::vector<std::vector<RealGradient>>>(
85 : "dPorousFlow_grad_porepressure_qp_dvar")),
86 2602 : _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
87 2602 : _num_phases(_dictator.numPhases()),
88 5204 : _gravity(this->template getParam<RealVectorValue>("gravity")),
89 2602 : _perm_derivs(_dictator.usePermDerivs()),
90 5204 : _full_upwind_threshold(this->template getParam<unsigned>("full_upwind_threshold")),
91 2602 : _fallback_scheme(
92 2602 : this->template getParam<MooseEnum>("fallback_scheme").template getEnum<FallbackEnum>()),
93 2602 : _proto_flux(_num_phases),
94 2602 : _jacobian(_num_phases),
95 2602 : _num_upwinds(),
96 2602 : _num_downwinds()
97 : {
98 : // The full-upwinding scheme identifies each element test function with a mesh node, which is
99 : // only valid for nodal (Lagrange) variables.
100 2602 : if (!_var.isNodal())
101 2 : mooseError("The variable '",
102 2 : _var.name(),
103 : "' is not a nodal (Lagrange) variable. This kernel uses full upwinding, which "
104 : "requires a nodal variable. For non-nodal variables use the non-upwinded "
105 : "PorousFlowFullySaturated* kernels or Kuzmin-Turek (KT) stabilisation instead.");
106 :
107 : #ifdef LIBMESH_HAVE_TBB_API
108 : if (libMesh::n_threads() > 1)
109 : mooseWarning("PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed "
110 : "accurately when using TBB and greater than 1 thread");
111 : #endif
112 2600 : }
113 :
114 : template <bool is_ad>
115 : void
116 25002 : PorousFlowDarcyBaseTempl<is_ad>::timestepSetup()
117 : {
118 25002 : GenericKernel<is_ad>::timestepSetup();
119 25002 : _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
120 25002 : _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
121 25002 : }
122 :
123 : template <bool is_ad>
124 : void
125 58270 : PorousFlowDarcyBaseTempl<is_ad>::jacobianSetup()
126 : {
127 58270 : GenericKernel<is_ad>::jacobianSetup();
128 58270 : _my_elem_darcy = nullptr;
129 58270 : }
130 :
131 : template <bool is_ad>
132 : GenericReal<is_ad>
133 199252576 : PorousFlowDarcyBaseTempl<is_ad>::darcyQp(unsigned int ph) const
134 : {
135 199252576 : return _grad_test[_i][_qp] *
136 199252576 : (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity));
137 : }
138 :
139 : template <bool is_ad>
140 : Real
141 648261696 : PorousFlowDarcyBaseTempl<is_ad>::darcyQpJacobian(unsigned int jvar, unsigned int ph) const
142 : {
143 : if constexpr (!is_ad)
144 : {
145 648261696 : if (_dictator.notPorousFlowVariable(jvar))
146 : return 0.0;
147 :
148 648261696 : const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
149 :
150 : RealVectorValue deriv =
151 648261696 : _permeability[_qp] * (_grad_phi[_j][_qp] * (*_dgrad_p_dgrad_var)[_qp][ph][pvar] -
152 648261696 : _phi[_j][_qp] * (*_dfluid_density_qp_dvar)[_qp][ph][pvar] * _gravity);
153 :
154 648261696 : deriv += _permeability[_qp] * ((*_dgrad_p_dvar)[_qp][ph][pvar] * _phi[_j][_qp]);
155 :
156 648261696 : if (_perm_derivs)
157 : {
158 4853504 : deriv += (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp] *
159 4853504 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
160 19414016 : for (const auto i : make_range(Moose::dim))
161 14560512 : deriv += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
162 14560512 : (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
163 : }
164 :
165 648261696 : return _grad_test[_i][_qp] * deriv;
166 : }
167 : else
168 : libmesh_ignore(jvar, ph);
169 0 : return 0.0;
170 : }
171 :
172 : template <bool is_ad>
173 : GenericReal<is_ad>
174 0 : PorousFlowDarcyBaseTempl<is_ad>::computeQpResidual()
175 : {
176 0 : mooseError("PorousFlowDarcyBase: computeQpResidual called");
177 : return 0.0;
178 : }
179 :
180 : template <bool is_ad>
181 : void
182 3542878 : PorousFlowDarcyBaseTempl<is_ad>::computeResidual()
183 : {
184 : if constexpr (!is_ad)
185 3542302 : computeResidualAndJacobian(JacRes::CALCULATE_RESIDUAL, 0);
186 : else
187 : {
188 576 : adComputeProtoFlux(false);
189 : // Assemble the real-valued residual from the value parts of the ADReal proto fluxes
190 576 : assembleProtoFluxResidual();
191 : }
192 3542878 : }
193 :
194 : template <bool is_ad>
195 : void
196 3542878 : PorousFlowDarcyBaseTempl<is_ad>::assembleProtoFluxResidual()
197 : {
198 3542878 : this->prepareVectorTag(this->_assembly, _var.number());
199 16481322 : for (_i = 0; _i < _test.size(); _i++)
200 27818872 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
201 14880428 : this->_local_re(_i) += MetaPhysicL::raw_value(_proto_flux[ph][_i]);
202 3542878 : this->accumulateTaggedLocalResidual();
203 :
204 3542878 : if (this->_has_save_in)
205 0 : for (unsigned int i = 0; i < this->_save_in.size(); i++)
206 0 : this->_save_in[i]->sys().solution().add_vector(this->_local_re,
207 0 : this->_save_in[i]->dofIndices());
208 3542878 : }
209 :
210 : template <bool is_ad>
211 : void
212 336 : PorousFlowDarcyBaseTempl<is_ad>::computeJacobian()
213 : {
214 : if constexpr (!is_ad)
215 336 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, _var.number());
216 : else
217 : // Block-diagonal (e.g. PJFNK) assembly path: this is the only Jacobian call we get.
218 0 : adComputeJacobian();
219 336 : }
220 :
221 : template <bool is_ad>
222 : void
223 48 : PorousFlowDarcyBaseTempl<is_ad>::adComputeJacobian()
224 : {
225 : if constexpr (is_ad)
226 : {
227 : // Compute ADReal proto fluxes with upwind counting, then extract the Jacobian.
228 : // Every block (diagonal and off-diagonal) is encoded in the ADReal derivatives, so a
229 : // single call to addJacobianWithoutConstraints captures all variable sensitivities.
230 : //
231 : // See PorousFlowLumpedKernelBase.C::computeJacobian for the detailed rationale of why
232 : // addJacobianWithoutConstraints (rather than the default ADKernel addJacobian) is required:
233 : // each residual row's column set must be read from its own derivatives instead of being shared
234 : // from row 0. There it is forced by mass-lumped nodal materials; here it is forced by the
235 : // upwinding scheme, which makes each node's proto flux depend on a different set of nodal DOFs.
236 48 : adComputeProtoFlux(true);
237 :
238 48 : const unsigned int num_nodes = _test.size();
239 48 : std::vector<ADReal> darcy_residuals(num_nodes, 0.0);
240 192 : for (const auto n : make_range(num_nodes))
241 288 : for (const auto ph : make_range(_num_phases))
242 144 : darcy_residuals[n] += _proto_flux[ph][n];
243 :
244 48 : this->addJacobianWithoutConstraints(
245 48 : this->_assembly, darcy_residuals, this->dofIndices(), _var.scalingFactor());
246 48 : }
247 48 : }
248 :
249 : template <bool is_ad>
250 : void
251 4281414 : PorousFlowDarcyBaseTempl<is_ad>::computeOffDiagJacobian(const unsigned int jvar)
252 : {
253 : if constexpr (!is_ad)
254 4281302 : computeResidualAndJacobian(JacRes::CALCULATE_JACOBIAN, jvar);
255 : else
256 : {
257 : libmesh_ignore(jvar);
258 : // Full (SMP) assembly calls this once per coupled variable; the first call performs the
259 : // single AD pass that fills every block, and the guard makes the rest no-ops. This mirrors
260 : // ADKernel::computeOffDiagJacobian and avoids the ADD-semantics double-counting.
261 112 : if (_my_elem_darcy != this->_current_elem)
262 : {
263 48 : adComputeJacobian();
264 48 : _my_elem_darcy = this->_current_elem;
265 : }
266 : }
267 4281414 : }
268 :
269 : template <bool is_ad>
270 : void
271 7824564 : PorousFlowDarcyBaseTempl<is_ad>::computeProtoFluxWithoutMobility()
272 : {
273 7824564 : const unsigned int num_nodes = _test.size();
274 17403721 : for (unsigned ph = 0; ph < _num_phases; ++ph)
275 : {
276 9579157 : _proto_flux[ph].assign(num_nodes, 0.0);
277 46040581 : for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
278 : {
279 36461424 : const Real jxw_coord = this->_JxW[_qp] * this->_coord[_qp];
280 235714000 : for (_i = 0; _i < num_nodes; ++_i)
281 199259392 : _proto_flux[ph][_i] += jxw_coord * darcyQp(ph);
282 : }
283 : }
284 7824564 : }
285 :
286 : template <bool is_ad>
287 : void
288 7824564 : PorousFlowDarcyBaseTempl<is_ad>::initializeUpwindTracking(unsigned elem, unsigned int num_nodes)
289 : {
290 7824564 : if (_num_upwinds.find(elem) == _num_upwinds.end())
291 : {
292 1793334 : _num_upwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
293 1793334 : _num_downwinds[elem] = std::vector<std::vector<unsigned>>(_num_phases);
294 1916185 : for (unsigned ph = 0; ph < _num_phases; ++ph)
295 : {
296 1019518 : _num_upwinds[elem][ph].assign(num_nodes, 0);
297 1019518 : _num_downwinds[elem][ph].assign(num_nodes, 0);
298 : }
299 : }
300 7824564 : }
301 :
302 : template <bool is_ad>
303 : void
304 2073075 : PorousFlowDarcyBaseTempl<is_ad>::updateUpwindCounts(unsigned elem, unsigned int num_nodes)
305 : {
306 4606532 : for (unsigned ph = 0; ph < _num_phases; ++ph)
307 11391251 : for (unsigned nod = 0; nod < num_nodes; ++nod)
308 : {
309 8857794 : if (_proto_flux[ph][nod] > 0)
310 4291519 : _num_upwinds[elem][ph][nod]++;
311 4566275 : else if (_proto_flux[ph][nod] < 0)
312 4235731 : _num_downwinds[elem][ph][nod]++;
313 : }
314 2073075 : }
315 :
316 : template <bool is_ad>
317 : std::vector<unsigned>
318 7824564 : PorousFlowDarcyBaseTempl<is_ad>::computeMaxSwaps(unsigned elem, unsigned int num_nodes) const
319 : {
320 7824564 : std::vector<unsigned> max_swaps(_num_phases, 0);
321 17403721 : for (unsigned ph = 0; ph < _num_phases; ++ph)
322 43097329 : for (unsigned nod = 0; nod < num_nodes; ++nod)
323 33518172 : max_swaps[ph] =
324 : std::max(max_swaps[ph],
325 33518172 : std::min(_num_upwinds.at(elem)[ph][nod], _num_downwinds.at(elem)[ph][nod]));
326 7824564 : return max_swaps;
327 0 : }
328 :
329 : template <bool is_ad>
330 : void
331 7824564 : PorousFlowDarcyBaseTempl<is_ad>::applyUpwinding(const std::vector<unsigned> & max_swaps,
332 : JacRes res_or_jac,
333 : unsigned int pvar)
334 : {
335 17403721 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
336 : {
337 9579157 : if (max_swaps[ph] < _full_upwind_threshold)
338 9577925 : fullyUpwind(res_or_jac, ph, pvar);
339 : else
340 : {
341 1232 : switch (_fallback_scheme)
342 : {
343 432 : case FallbackEnum::QUICK:
344 432 : quickUpwind(res_or_jac, ph, pvar);
345 432 : break;
346 800 : case FallbackEnum::HARMONIC:
347 800 : harmonicMean(res_or_jac, ph, pvar);
348 800 : break;
349 : }
350 : }
351 : }
352 7824564 : }
353 :
354 : template <bool is_ad>
355 : void
356 624 : PorousFlowDarcyBaseTempl<is_ad>::adComputeProtoFlux(bool do_counting)
357 : {
358 624 : const unsigned int num_nodes = _test.size();
359 624 : computeProtoFluxWithoutMobility();
360 :
361 : // Initialise upwind-tracking maps for this element on first encounter
362 624 : const unsigned elem = this->_current_elem->id();
363 624 : initializeUpwindTracking(elem, num_nodes);
364 :
365 624 : if (do_counting)
366 48 : updateUpwindCounts(elem, num_nodes);
367 :
368 : // Determine how many upwind-downwind swaps have occurred this timestep
369 624 : const std::vector<unsigned> max_swaps = computeMaxSwaps(elem, num_nodes);
370 :
371 : // Apply mobility via the chosen upwinding scheme (always in residual mode for AD)
372 624 : applyUpwinding(max_swaps, JacRes::CALCULATE_RESIDUAL, 0);
373 624 : }
374 :
375 : template <bool is_ad>
376 : void
377 7823940 : PorousFlowDarcyBaseTempl<is_ad>::computeResidualAndJacobian(JacRes res_or_jac, unsigned int jvar)
378 : {
379 7823940 : if ((res_or_jac == JacRes::CALCULATE_JACOBIAN) && _dictator.notPorousFlowVariable(jvar))
380 0 : return;
381 :
382 : // The PorousFlow variable index corresponding to the variable number jvar
383 : const unsigned int pvar =
384 7823940 : ((res_or_jac == JacRes::CALCULATE_JACOBIAN) ? _dictator.porousFlowVariableNum(jvar) : 0);
385 :
386 7823940 : this->prepareMatrixTag(this->_assembly, _var.number(), jvar);
387 7823940 : if ((this->_local_ke.n() == 0) &&
388 : (res_or_jac == JacRes::CALCULATE_JACOBIAN)) // this removes a problem
389 : // encountered in the
390 : // initial timestep when
391 : // use_displaced_mesh=true
392 : return;
393 :
394 : // The number of nodes in the element
395 7823940 : const unsigned int num_nodes = _test.size();
396 :
397 : // Compute the residual and jacobian without the mobility terms. Even if we are computing the
398 : // Jacobian we still need this in order to see which nodes are upwind and which are downwind.
399 7823940 : computeProtoFluxWithoutMobility();
400 :
401 : // for this element, record whether each node is "upwind" or "downwind" (or neither)
402 7823940 : const unsigned elem = this->_current_elem->id();
403 7823940 : initializeUpwindTracking(elem, num_nodes);
404 : // record the information once per nonlinear iteration
405 7823940 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN && jvar == _var.number())
406 2073027 : updateUpwindCounts(elem, num_nodes);
407 :
408 : // based on _num_upwinds and _num_downwinds, calculate the maximum number
409 : // of upwind-downwind swaps that have been encountered in this timestep
410 : // for this element
411 7823940 : const std::vector<unsigned> max_swaps = computeMaxSwaps(elem, num_nodes);
412 :
413 : // size the _jacobian correctly and calculate it for the case residual = _proto_flux
414 7823940 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
415 : {
416 9634140 : for (unsigned ph = 0; ph < _num_phases; ++ph)
417 : {
418 5352502 : _jacobian[ph].resize(this->_local_ke.m());
419 23990102 : for (_i = 0; _i < _test.size(); _i++)
420 : {
421 18637600 : _jacobian[ph][_i].assign(this->_local_ke.n(), 0.0);
422 101594820 : for (_j = 0; _j < _phi.size(); _j++)
423 731218916 : for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
424 648261696 : _jacobian[ph][_i][_j] +=
425 648261696 : this->_JxW[_qp] * this->_coord[_qp] * darcyQpJacobian(jvar, ph);
426 : }
427 : }
428 : }
429 :
430 : // Loop over all the phases, computing the mass flux, which
431 : // gets placed into _proto_flux, and the derivative of this
432 : // which gets placed into _jacobian
433 7823940 : applyUpwinding(max_swaps, res_or_jac, pvar);
434 :
435 : // Add results to the Residual or Jacobian
436 7823940 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
437 3542302 : assembleProtoFluxResidual();
438 :
439 7823940 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
440 : {
441 20040182 : for (_i = 0; _i < _test.size(); _i++)
442 90004884 : for (_j = 0; _j < _phi.size(); _j++)
443 157203560 : for (unsigned int ph = 0; ph < _num_phases; ++ph)
444 82957220 : this->_local_ke(_i, _j) += _jacobian[ph][_i][_j];
445 :
446 4281638 : this->accumulateTaggedLocalMatrix();
447 :
448 4281638 : if (this->_has_diag_save_in && jvar == _var.number())
449 : {
450 : unsigned int rows = this->_local_ke.m();
451 0 : DenseVector<Number> diag(rows);
452 0 : for (unsigned int i = 0; i < rows; i++)
453 0 : diag(i) = this->_local_ke(i, i);
454 :
455 0 : for (unsigned int i = 0; i < this->_diag_save_in.size(); i++)
456 0 : this->_diag_save_in[i]->sys().solution().add_vector(diag,
457 0 : this->_diag_save_in[i]->dofIndices());
458 : }
459 : }
460 7823940 : }
461 :
462 : template <bool is_ad>
463 : void
464 9577925 : PorousFlowDarcyBaseTempl<is_ad>::fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
465 : {
466 : // res_or_jac and pvar drive the hand-coded non-AD Jacobian only; the AD path ignores them
467 : if constexpr (is_ad)
468 : libmesh_ignore(res_or_jac, pvar);
469 :
470 : /**
471 : * Perform the full upwinding by multiplying the residuals at the upstream nodes by their
472 : * mobilities.
473 : * Mobility is different for each phase, and in each situation:
474 : * mobility = density / viscosity for single-component Darcy flow
475 : * mobility = mass_fraction * density * relative_perm / viscosity for multi-component,
476 : *multiphase flow
477 : * mobility = enthalpy * density * relative_perm / viscosity for heat convection
478 : *
479 : * The residual for the kernel is the sum over Darcy fluxes for each phase.
480 : * The Darcy flux for a particular phase is
481 : * R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
482 : * for node i. where int is the integral over the element.
483 : * However, in fully-upwind, the first step is to take the mobility outside the integral,
484 : * which was done in the _proto_flux calculation above.
485 : *
486 : * NOTE: Physically _proto_flux[i][ph] is a measure of fluid of phase ph flowing out of node i.
487 : * If we had left in mobility, it would be exactly the component mass flux flowing out of node
488 : *i.
489 : *
490 : * This leads to the definition of upwinding:
491 : *
492 : * If _proto_flux(i)[ph] is positive then we use mobility_i. That is we use the upwind value of
493 : * mobility.
494 : *
495 : * The final subtle thing is we must also conserve fluid mass: the total component mass flowing
496 : * out of node i must be the sum of the masses flowing into the other nodes.
497 : **/
498 :
499 : // The number of nodes in the element
500 9577925 : const unsigned int num_nodes = _test.size();
501 :
502 : GenericReal<is_ad> mob;
503 : // Define variables used to ensure mass conservation
504 624 : GenericReal<is_ad> total_mass_out = 0.0;
505 624 : GenericReal<is_ad> total_in = 0.0;
506 :
507 : // The following holds derivatives of these (non-AD path only)
508 : std::vector<Real> dtotal_mass_out;
509 : std::vector<Real> dtotal_in;
510 : if constexpr (!is_ad)
511 9577301 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
512 : {
513 5352342 : dtotal_mass_out.assign(num_nodes, 0.0);
514 5352342 : dtotal_in.assign(num_nodes, 0.0);
515 : }
516 :
517 : // Perform the upwinding using the mobility
518 9577925 : std::vector<bool> upwind_node(num_nodes);
519 43091233 : for (unsigned int n = 0; n < num_nodes; ++n)
520 : {
521 33513308 : if (_proto_flux[ph][n] >= 0.0) // upstream node
522 : {
523 : upwind_node[n] = true;
524 : // The mobility at the upstream node
525 17371159 : mob = mobility(n, ph);
526 : if constexpr (!is_ad)
527 17370179 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
528 : {
529 : // The derivative of the mobility wrt the PorousFlow variable
530 9698873 : const Real dmob = dmobility(n, ph, pvar);
531 :
532 53005013 : for (_j = 0; _j < _phi.size(); _j++)
533 43306140 : _jacobian[ph][n][_j] *= mob;
534 :
535 9698873 : if (_test.size() == _phi.size())
536 : /* mobility at node=n depends only on the variables at node=n, by construction. For
537 : * linear-lagrange variables, this means that Jacobian entries involving the derivative
538 : * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
539 : * [n][n] in the line below. However, for other variable types (eg constant monomials)
540 : * I cannot tell what variable number contributes to the derivative. However, in all
541 : * cases I can possibly imagine, the derivative is zero anyway, since in the full
542 : * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
543 : */
544 8458326 : _jacobian[ph][n][n] += dmob * MetaPhysicL::raw_value(_proto_flux[ph][n]);
545 :
546 53005013 : for (_j = 0; _j < _phi.size(); _j++)
547 43306140 : dtotal_mass_out[_j] += _jacobian[ph][n][_j];
548 : }
549 17371159 : _proto_flux[ph][n] *= mob;
550 17371159 : total_mass_out += _proto_flux[ph][n];
551 : }
552 : else
553 : {
554 : upwind_node[n] = false;
555 16142149 : total_in -= _proto_flux[ph][n]; /// note the -= means the result is positive
556 : if constexpr (!is_ad)
557 16141161 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
558 48586831 : for (_j = 0; _j < _phi.size(); _j++)
559 39648712 : dtotal_in[_j] -= _jacobian[ph][n][_j];
560 : }
561 : }
562 :
563 : // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
564 : // weighted by their proto_flux values
565 43091233 : for (unsigned int n = 0; n < num_nodes; ++n)
566 : {
567 33513308 : if (!upwind_node[n]) // downstream node
568 : {
569 : if constexpr (!is_ad)
570 16141161 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
571 48586831 : for (_j = 0; _j < _phi.size(); _j++)
572 : {
573 39648712 : _jacobian[ph][n][_j] *= MetaPhysicL::raw_value(total_mass_out / total_in);
574 39648712 : _jacobian[ph][n][_j] +=
575 39648712 : MetaPhysicL::raw_value(_proto_flux[ph][n]) *
576 39648712 : (dtotal_mass_out[_j] / MetaPhysicL::raw_value(total_in) -
577 39648712 : dtotal_in[_j] * MetaPhysicL::raw_value(total_mass_out) /
578 39648712 : MetaPhysicL::raw_value(total_in) / MetaPhysicL::raw_value(total_in));
579 : }
580 16142149 : _proto_flux[ph][n] *= total_mass_out / total_in;
581 : }
582 : }
583 9577925 : }
584 :
585 : template <bool is_ad>
586 : void
587 432 : PorousFlowDarcyBaseTempl<is_ad>::quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
588 : {
589 : // res_or_jac and pvar drive the hand-coded non-AD Jacobian only; the AD path ignores them
590 : if constexpr (is_ad)
591 : libmesh_ignore(res_or_jac, pvar);
592 :
593 : // The number of nodes in the element
594 432 : const unsigned int num_nodes = _test.size();
595 :
596 : // Use the raw nodal mobility
597 2096 : for (unsigned int n = 0; n < num_nodes; ++n)
598 : {
599 : // The mobility at the node
600 1664 : const GenericReal<is_ad> mob = mobility(n, ph);
601 : if constexpr (!is_ad)
602 1664 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
603 : {
604 : // The derivative of the mobility wrt the PorousFlow variable
605 224 : const Real dmob = dmobility(n, ph, pvar);
606 :
607 1056 : for (_j = 0; _j < _phi.size(); _j++)
608 832 : _jacobian[ph][n][_j] *= mob;
609 :
610 224 : if (_test.size() == _phi.size())
611 : /* mobility at node=n depends only on the variables at node=n, by construction. For
612 : * linear-lagrange variables, this means that Jacobian entries involving the derivative
613 : * of mobility will only be nonzero for derivatives wrt variables at node=n. Hence the
614 : * [n][n] in the line below. However, for other variable types (eg constant monomials)
615 : * I cannot tell what variable number contributes to the derivative. However, in all
616 : * cases I can possibly imagine, the derivative is zero anyway, since in the full
617 : * upwinding scheme, mobility shouldn't depend on these other sorts of variables.
618 : */
619 224 : _jacobian[ph][n][n] += dmob * MetaPhysicL::raw_value(_proto_flux[ph][n]);
620 : }
621 1664 : _proto_flux[ph][n] *= mob;
622 : }
623 432 : }
624 :
625 : template <bool is_ad>
626 : void
627 800 : PorousFlowDarcyBaseTempl<is_ad>::harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
628 : {
629 : // res_or_jac and pvar drive the hand-coded non-AD Jacobian only; the AD path ignores them
630 : if constexpr (is_ad)
631 : libmesh_ignore(res_or_jac, pvar);
632 :
633 : // The number of nodes in the element
634 800 : const unsigned int num_nodes = _test.size();
635 :
636 800 : std::vector<GenericReal<is_ad>> mob(num_nodes);
637 : unsigned num_zero = 0;
638 0 : GenericReal<is_ad> harmonic_mob = 0;
639 4000 : for (unsigned n = 0; n < num_nodes; ++n)
640 : {
641 3200 : mob[n] = mobility(n, ph);
642 3200 : if (MetaPhysicL::raw_value(mob[n]) == 0.0)
643 : {
644 0 : num_zero++;
645 : }
646 : else
647 3200 : harmonic_mob += 1.0 / mob[n];
648 : }
649 800 : if (num_zero > 0)
650 0 : harmonic_mob = 0.0;
651 : else
652 800 : harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
653 :
654 : // d(harmonic_mob)/d(PorousFlow variable at node n) -- non-AD path only
655 800 : std::vector<Real> dharmonic_mob(num_nodes, 0.0);
656 : if constexpr (!is_ad)
657 800 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
658 : {
659 96 : const Real harm2 = std::pow(MetaPhysicL::raw_value(harmonic_mob), 2) / (1.0 * num_nodes);
660 96 : if (num_zero == 0)
661 480 : for (unsigned n = 0; n < num_nodes; ++n)
662 384 : dharmonic_mob[n] =
663 384 : dmobility(n, ph, pvar) * harm2 / std::pow(MetaPhysicL::raw_value(mob[n]), 2);
664 0 : else if (num_zero == 1)
665 0 : for (unsigned n = 0; n < num_nodes; ++n)
666 0 : if (MetaPhysicL::raw_value(mob[n]) == 0.0)
667 : {
668 0 : dharmonic_mob[n] = num_nodes * dmobility(n, ph, pvar); // other derivs are zero
669 0 : break;
670 : }
671 : // if num_zero > 1 then all dharmonic_mob = 0.0
672 : }
673 :
674 : if constexpr (!is_ad)
675 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
676 480 : for (unsigned n = 0; n < num_nodes; ++n)
677 1920 : for (_j = 0; _j < _phi.size(); _j++)
678 : {
679 1536 : _jacobian[ph][n][_j] *= MetaPhysicL::raw_value(harmonic_mob);
680 1536 : if (_test.size() == _phi.size())
681 1536 : _jacobian[ph][n][_j] += dharmonic_mob[_j] * MetaPhysicL::raw_value(_proto_flux[ph][n]);
682 : }
683 :
684 4000 : for (unsigned n = 0; n < num_nodes; ++n)
685 3200 : _proto_flux[ph][n] *= harmonic_mob;
686 800 : }
687 :
688 : template <bool is_ad>
689 : GenericReal<is_ad>
690 0 : PorousFlowDarcyBaseTempl<is_ad>::mobility(unsigned nodenum, unsigned phase) const
691 : {
692 0 : return _fluid_density_node[nodenum][phase] / _fluid_viscosity[nodenum][phase];
693 : }
694 :
695 : template <bool is_ad>
696 : Real
697 0 : PorousFlowDarcyBaseTempl<is_ad>::dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
698 : {
699 : if constexpr (!is_ad)
700 : {
701 0 : Real dm = (*_dfluid_density_node_dvar)[nodenum][phase][pvar] / _fluid_viscosity[nodenum][phase];
702 0 : dm -= _fluid_density_node[nodenum][phase] * (*_dfluid_viscosity_dvar)[nodenum][phase][pvar] /
703 : std::pow(_fluid_viscosity[nodenum][phase], 2);
704 : return dm;
705 : }
706 : else
707 : libmesh_ignore(nodenum, phase, pvar);
708 0 : return 0.0;
709 : }
710 :
711 : template class PorousFlowDarcyBaseTempl<false>;
712 : template class PorousFlowDarcyBaseTempl<true>;
|
