17 #include "libmesh/quadrature.h" 24 "Gravitational acceleration vector downwards (m/s^2)");
26 "PorousFlowDictator",
"The UserObject that holds the list of PorousFlow variable names");
27 params.
addParam<
unsigned>(
"full_upwind_threshold",
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 MooseEnum fallback_enum(
"quick harmonic",
"quick");
36 "quick: use nodal mobility without " 37 "preserving mass. harmonic: use a " 38 "harmonic mean of nodal mobilities " 39 "and preserve fluid mass");
46 _permeability(getMaterialProperty<
RealTensorValue>(
"PorousFlow_permeability_qp")),
48 getMaterialProperty<
std::vector<
RealTensorValue>>(
"dPorousFlow_permeability_qp_dvar")),
50 "dPorousFlow_permeability_qp_dgradvar")),
52 getMaterialProperty<
std::vector<
Real>>(
"PorousFlow_fluid_phase_density_nodal")),
53 _dfluid_density_node_dvar(getMaterialProperty<
std::vector<
std::vector<
Real>>>(
54 "dPorousFlow_fluid_phase_density_nodal_dvar")),
55 _fluid_density_qp(getMaterialProperty<
std::vector<
Real>>(
"PorousFlow_fluid_phase_density_qp")),
56 _dfluid_density_qp_dvar(getMaterialProperty<
std::vector<
std::vector<
Real>>>(
57 "dPorousFlow_fluid_phase_density_qp_dvar")),
58 _fluid_viscosity(getMaterialProperty<
std::vector<
Real>>(
"PorousFlow_viscosity_nodal")),
59 _dfluid_viscosity_dvar(
60 getMaterialProperty<
std::vector<
std::vector<
Real>>>(
"dPorousFlow_viscosity_nodal_dvar")),
61 _pp(getMaterialProperty<
std::vector<
Real>>(
"PorousFlow_porepressure_nodal")),
62 _grad_p(getMaterialProperty<
std::vector<
RealGradient>>(
"PorousFlow_grad_porepressure_qp")),
63 _dgrad_p_dgrad_var(getMaterialProperty<
std::vector<
std::vector<
Real>>>(
64 "dPorousFlow_grad_porepressure_qp_dgradvar")),
66 "dPorousFlow_grad_porepressure_qp_dvar")),
68 _num_phases(_dictator.numPhases()),
70 _perm_derivs(_dictator.usePermDerivs()),
71 _full_upwind_threshold(getParam<unsigned>(
"full_upwind_threshold")),
73 _proto_flux(_num_phases),
74 _jacobian(_num_phases),
78 #ifdef LIBMESH_HAVE_TBB_API 80 mooseWarning(
"PorousFlowDarcyBase: num_upwinds and num_downwinds may not be computed " 81 "accurately when using TBB and greater than 1 thread");
89 _num_upwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
90 _num_downwinds = std::unordered_map<unsigned, std::vector<std::vector<unsigned>>>();
129 mooseError(
"PorousFlowDarcyBase: computeQpResidual called");
158 const unsigned int pvar =
169 const unsigned int num_nodes =
_test.size();
178 for (
_i = 0;
_i < num_nodes; ++
_i)
200 for (
unsigned nod = 0; nod < num_nodes; ++nod)
216 for (
unsigned nod = 0; nod < num_nodes; ++nod)
217 max_swaps[ph] = std::max(
268 for (
unsigned int i = 0; i <
_save_in.size(); i++)
285 for (
unsigned int i = 0; i < rows; i++)
327 const unsigned int num_nodes =
_test.size();
332 Real total_mass_out = 0.0;
336 std::vector<Real> dtotal_mass_out;
337 std::vector<Real> dtotal_in;
340 dtotal_mass_out.assign(num_nodes, 0.0);
341 dtotal_in.assign(num_nodes, 0.0);
345 std::vector<bool> upwind_node(num_nodes);
346 for (
unsigned int n = 0; n < num_nodes; ++n)
350 upwind_node[n] =
true;
380 upwind_node[n] =
false;
390 for (
unsigned int n = 0; n < num_nodes; ++n)
400 dtotal_in[
_j] * total_mass_out / total_in / total_in);
411 const unsigned int num_nodes =
_test.size();
417 for (
unsigned int n = 0; n < num_nodes; ++n)
448 const unsigned int num_nodes =
_test.size();
450 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 for (
unsigned n = 0; n < num_nodes; ++n)
459 zero_mobility_node = n;
463 harmonic_mob += 1.0 / mob[n];
468 harmonic_mob = (1.0 * num_nodes) / harmonic_mob;
471 std::vector<Real> dharmonic_mob(num_nodes, 0.0);
474 const Real harm2 =
std::pow(harmonic_mob, 2) / (1.0 * num_nodes);
476 for (
unsigned n = 0; n < num_nodes; ++n)
478 else if (num_zero == 1)
479 dharmonic_mob[zero_mobility_node] =
480 num_nodes *
dmobility(zero_mobility_node, ph, pvar);
485 for (
unsigned n = 0; n < num_nodes; ++n)
494 for (
unsigned n = 0; n < num_nodes; ++n)
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_node_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables (at the node) ...
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dpermeability_dgradvar
d(permeabiity)/d(grad(PorousFlow variable))
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_qp_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables (at the qp) ...
const MaterialProperty< std::vector< RealTensorValue > > & _dpermeability_dvar
d(permeabiity)/d(PorousFlow variable)
static InputParameters validParams()
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
void accumulateTaggedLocalResidual()
std::vector< MooseVariableFEBase *> _save_in
const MooseArray< Real > & _JxW
enum PorousFlowDarcyBase::FallbackEnum _fallback_scheme
unsigned int number() const
const VariablePhiGradient & _grad_phi
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< RealTensorValue > & _permeability
Permeability of porous material.
FallbackEnum
If full upwinding is failing due to nodes swapping between upwind and downwind in successive nonlinea...
std::unordered_map< unsigned, std::vector< std::vector< unsigned > > > _num_downwinds
Number of nonlinear iterations (in this timestep and this element) that a node is an downwind node fo...
std::vector< std::vector< Real > > _proto_flux
The Darcy flux.
static constexpr std::size_t dim
virtual void computeResidualAndJacobian() override
std::vector< std::vector< std::vector< Real > > > _jacobian
Derivative of _proto_flux with respect to nodal variables.
virtual Real darcyQpJacobian(unsigned int jvar, unsigned int ph) const
Jacobian of the Darcy part of the flux.
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
DenseMatrix< Number > _local_ke
void mooseWarning(Args &&... args) const
virtual void computeResidual() override
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
TensorValue< Real > RealTensorValue
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
Derivative of the fluid viscosity for each phase wrt PorousFlow variables.
const VariableTestValue & _test
virtual void timestepSetup()
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
static InputParameters validParams()
std::vector< MooseVariableFEBase *> _diag_save_in
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
const unsigned int _num_phases
The number of fluid phases.
const QBase *const & _qrule
virtual void computeOffDiagJacobian(unsigned int jvar) override
std::unordered_map< unsigned, std::vector< std::vector< unsigned > > > _num_upwinds
Number of nonlinear iterations (in this timestep and this element) that a node is an upwind node for ...
void accumulateTaggedLocalMatrix()
void harmonicMean(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian by using the harmonic mean of the nodal mobilities for the entire ...
virtual Real computeQpResidual() override
virtual void timestepSetup() override
const MooseArray< Real > & _coord
const RealVectorValue _gravity
Gravity. Defaults to 9.81 m/s^2.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real mobility(unsigned nodenum, unsigned phase) const
The mobility of the fluid.
PorousFlowDarcyBase(const InputParameters ¶meters)
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)
void quickUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using the nodal mobilities, but without conserving fluid mass...
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const VariableTestGradient & _grad_test
virtual Real dmobility(unsigned nodenum, unsigned phase, unsigned pvar) const
The derivative of mobility with respect to PorousFlow variable pvar.
const unsigned _full_upwind_threshold
If the number of upwind-downwind swaps is less than this amount then full upwinding is used...
DenseVector< Number > _local_re
virtual void computeJacobian() override
IntRange< T > make_range(T beg, T end)
void mooseError(Args &&... args) const
const Elem *const & _current_elem
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
virtual Real darcyQp(unsigned int ph) const
The Darcy part of the flux (this is the non-upwinded part)
const VariablePhiValue & _phi
const bool _perm_derivs
Flag to check whether permeabiity derivatives are non-zero.
MooseUnits pow(const MooseUnits &, int)
const MaterialProperty< std::vector< Real > > & _fluid_density_node
Fluid density for each phase (at the node)
void fullyUpwind(JacRes res_or_jac, unsigned int ph, unsigned int pvar)
Calculate the residual or Jacobian using full upwinding.