https://mooseframework.inl.gov
ConservativeAdvection.C
Go to the documentation of this file.
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 "ConservativeAdvection.h"
11 #include "SystemBase.h"
12 
13 using namespace libMesh;
14 
17 
18 template <bool is_ad>
21 {
23  params.addClassDescription("Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak "
24  "form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$. Velocity can be "
25  "given as 1) a variable, for which the gradient is automatically "
26  "taken, 2) a vector variable, or a 3) vector material.");
27  params.addParam<MaterialPropertyName>(
28  "velocity_scalar_coef",
29  "1.0",
30  "Name of material property multiplied against the velocity to scale advection strength.");
31  MooseEnum upwinding_type("none full", "none");
32  params.addParam<MooseEnum>("upwinding_type",
33  upwinding_type,
34  "Type of upwinding used. None: Typically results in overshoots and "
35  "undershoots, but numerical diffusion is minimized. Full: Overshoots "
36  "and undershoots are avoided, but numerical diffusion is large");
37  params.addParam<MaterialPropertyName>("advected_quantity",
38  "An optional material property to be advected. If not "
39  "supplied, then the variable will be used.");
40  params.addCoupledVar(
41  "velocity_as_variable_gradient",
42  "Gradient of this coupled variable is used to define the advection velocity. "
43  "Can be supplied instead of velocity material or velocity variable.");
44  return params;
45 }
46 
47 template <>
50 {
51  InputParameters params = generalParams();
52  params.addCoupledVar("velocity", "Velocity vector");
53  params.deprecateCoupledVar("velocity", "velocity_variable", "12/31/2025");
54  params.addParam<MaterialPropertyName>("velocity_material", "Velocity vector given as a material");
55  return params;
56 }
57 
58 template <>
61 {
62  InputParameters params = generalParams();
63  params.addCoupledVar("velocity_variable", "Velocity vector given as a variable");
64  params.addParam<MaterialPropertyName>("velocity", "Velocity vector given as a material");
65  params.deprecateParam("velocity", "velocity_material", "12/31/2025");
66  return params;
67 }
68 
69 template <bool is_ad>
71  : GenericKernel<is_ad>(parameters),
72  _scalar(this->template getGenericMaterialProperty<Real, is_ad>("velocity_scalar_coef")),
73  _coupled_variable_present(isParamValid("velocity_as_variable_gradient")),
74  _coupled_variable_var(_coupled_variable_present ? coupled("velocity_as_variable_gradient") : 0),
75  _velocity(
76  _coupled_variable_present
77  ? &this->template coupledGenericGradient<is_ad>("velocity_as_variable_gradient")
78  : (this->isParamValid("velocity_variable")
79  ? &this->template coupledGenericVectorValue<is_ad>("velocity_variable")
80  : (this->isParamValid("velocity_material")
81  ? &this->template getGenericMaterialProperty<RealVectorValue, is_ad>(
82  "velocity_material")
83  .get()
84  : nullptr))),
85  _user_supplied_adv_quant(isParamValid("advected_quantity")),
86  _adv_quant(
87  _user_supplied_adv_quant
88  ? this->template getGenericMaterialProperty<Real, is_ad>("advected_quantity").get()
89  : _u),
90  _upwinding(
91  this->template getParam<MooseEnum>("upwinding_type").template getEnum<UpwindingType>()),
92  _u_nodal(_var.template genericDofValues<is_ad>()),
93  _upwind_node(0),
94  _dtotal_mass_out(0)
95 {
97  paramError("velocity_as_variable_gradient",
98  "Use a different kernel (i.e., diffusion) if the gradient used as the velocity is "
99  "the same as the member variable");
100  if (_upwinding != UpwindingType::none && this->isParamValid("advected_quantity"))
101  paramError(
102  "advected_quantity",
103  "Upwinding is not compatible with an advected quantity that is not the primary variable.");
104 
105  if (!_velocity || (_coupled_variable_present && this->isParamValid("velocity_material")) ||
106  (_coupled_variable_present && this->isParamValid("velocity_variable")) ||
107  (this->isParamValid("velocity_variable") && this->isParamValid("velocity_material")))
108  paramError(
109  "velocity_as_variable_gradient",
110  "One and only one of the following input variables must be specified: velocity_variable, "
111  "velocity_material, or velocity_as_variable_gradient.");
112 
113  if (this->_has_diag_save_in)
114  paramError("diag_save_in",
115  "_local_ke not computed for global AD indexing. Save-in is deprecated anyway. Use "
116  "the tagging system instead.");
117 }
118 
119 template <bool is_ad>
122 {
123  return -_grad_test[_i][_qp] * (*_velocity)[_qp] * _scalar[_qp];
124 }
125 
126 template <bool is_ad>
129 {
130  // This is the no-upwinded version
131  // It gets called via GenericKernel<is_ad>::computeResidual()
132  return negSpeedQp() * _adv_quant[_qp];
133 }
134 
135 template <>
136 Real
138 {
139  // This is the no-upwinded version
140  // It gets called via GenericKernel<false>::computeJacobian()
141  if (!_user_supplied_adv_quant)
142  return negSpeedQp() * _phi[_j][_qp];
143  return 0.0;
144 }
145 
146 template <>
147 Real
149 {
150  mooseError("Internal error, should never get here when using AD");
151  return 0.0;
152 }
153 
154 template <>
155 Real
157 {
158  // This is the non-upwinded version
159  // It gets called via GenericKernel<false>::computeOffDiagJacobian()
160  if (_coupled_variable_present && _coupled_variable_var == jvar)
161  return -_grad_test[_i][_qp] * _grad_phi[_j][_qp] * _adv_quant[_qp] * _scalar[_qp];
162  else
163  return 0.0;
164 }
165 
166 template <>
167 Real
169 {
170  mooseError("Internal error, should never get here when using AD");
171  return 0.0;
172 }
173 
174 template <bool is_ad>
175 void
177 {
178  switch (_upwinding)
179  {
180  case UpwindingType::none:
182  break;
183  case UpwindingType::full:
184  fullUpwind(JacRes::CALCULATE_RESIDUAL);
185  break;
186  }
187 }
188 
189 template <bool is_ad>
190 void
192 {
193  switch (_upwinding)
194  {
195  case UpwindingType::none:
197  break;
198  case UpwindingType::full:
199  fullUpwind(JacRes::CALCULATE_JACOBIAN);
200  break;
201  }
202 }
203 
204 template <bool is_ad>
205 void
207 {
208  // The number of nodes in the element
209  const unsigned int num_nodes = _test.size();
210 
211  // Even if we are computing the Jacobian we still need to compute the outflow from each node to
212  // see which nodes are upwind and which are downwind
213  _my_local_re.resize(_var.dofIndices().size());
214 
215  if (!is_ad && (res_or_jac == JacRes::CALCULATE_JACOBIAN))
216  prepareMatrixTag(this->_assembly, _var.number(), _var.number());
217 
218  // Compute the outflux from each node and store in _my_local_re
219  // If _my_local_re is positive at the node, mass (or whatever the Variable represents) is flowing
220  // out of the node
221  _upwind_node.resize(num_nodes);
222  for (_i = 0; _i < num_nodes; ++_i)
223  {
224  for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
225  _my_local_re(_i) += this->_JxW[_qp] * this->_coord[_qp] * negSpeedQp();
226  _upwind_node[_i] = (MetaPhysicL::raw_value(_my_local_re(_i)) >= 0.0);
227  }
228 
229  // Variables used to ensure mass conservation
230  GenericReal<is_ad> total_mass_out = 0.0;
231  GenericReal<is_ad> total_in = 0.0;
232  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
233  _dtotal_mass_out.assign(num_nodes, 0.0);
234 
235  for (const auto n : make_range(num_nodes))
236  {
237  if (_upwind_node[n])
238  {
239  if constexpr (!is_ad)
240  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
241  {
242  if (_test.size() == _phi.size())
243  /* u at node=n depends only on the u at node=n, by construction. For
244  * linear-lagrange variables, this means that Jacobian entries involving the derivative
245  * will only be nonzero for derivatives wrt variable at node=n. Hence the
246  * (n, n) in the line below. The above "if" statement catches other variable types
247  * (eg constant monomials)
248  */
249  _local_ke(n, n) += _my_local_re(n);
250 
251  _dtotal_mass_out[n] += _local_ke(n, n);
252  }
253  _my_local_re(n) *= _u_nodal[n];
254  total_mass_out += _my_local_re(n);
255  }
256  else // downwind node
257  total_in -= _my_local_re(n); // note the -= means the result is positive
258  }
259 
260  // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
261  // weighted by their local_re values
262  for (const auto n : make_range(num_nodes))
263  if (!_upwind_node[n]) // downwind node
264  {
265  if constexpr (!is_ad)
266  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
267  for (_j = 0; _j < _phi.size(); _j++)
268  _local_ke(n, _j) += _my_local_re(n) * _dtotal_mass_out[_j] / total_in;
269  _my_local_re(n) *= total_mass_out / total_in;
270  }
271 
272  // Add the result to the residual and jacobian
273  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
274  {
275  this->addResiduals(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
276 
277  if (this->_has_save_in)
278  {
279  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
280  for (const auto & var : this->_save_in)
281  var->sys().solution().add_vector(MetaPhysicL::raw_value(_my_local_re), var->dofIndices());
282  }
283  }
284 
285  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
286  {
287  if constexpr (is_ad)
288  this->addJacobian(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
289  else
290  accumulateTaggedLocalMatrix();
291  }
292 }
293 
294 template class ConservativeAdvectionTempl<false>;
295 template class ConservativeAdvectionTempl<true>;
Moose::GenericType< Real, is_ad > GenericReal
Definition: MooseTypes.h:649
virtual void computeResidual() override
Compute this Kernel&#39;s contribution to the residual.
Definition: Kernel.C:92
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
For coupling standard variables.
virtual void computeResidual() override
Compute this Kernel&#39;s contribution to the residual.
const unsigned int _coupled_variable_var
Coupled variable variable number.
MooseVariable & _var
This is a regular kernel so we cast to a regular MooseVariable.
Definition: Kernel.h:72
void paramError(const std::string &param, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
Definition: MooseBase.h:439
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:323
unsigned int number() const
Get variable number coming from libMesh.
static InputParameters validParams()
Definition: GenericKernel.h:19
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
Definition: MooseUtils.h:1133
const MooseArray< GenericRealVectorValue< is_ad > > * _velocity
advection velocity
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
const bool _coupled_variable_present
Flag to determine if coupled variable is present.
enum ConservativeAdvectionTempl::UpwindingType _upwinding
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
bool _has_diag_save_in
The aux variables to save the diagonal Jacobian contributions to.
Definition: KernelBase.h:69
virtual GenericReal< is_ad > negSpeedQp() const
Returns - _grad_test * velocity.
registerMooseObject("MooseApp", ConservativeAdvection)
void deprecateParam(const std::string &old_name, const std::string &new_name, const std::string &removal_date)
static InputParameters generalParams()
virtual GenericReal< is_ad > computeQpResidual() override
Compute this Kernel&#39;s contribution to the residual at the current quadrature point.
ConservativeAdvectionTempl(const InputParameters &parameters)
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
Advection of the variable by the velocity provided by the user.
void deprecateCoupledVar(const std::string &old_name, const std::string &new_name, const std::string &removal_date)
void fullUpwind(JacRes res_or_jac)
Calculates the fully-upwind Residual and Jacobian (depending on res_or_jac)
virtual Real computeQpJacobian() override
Compute this Kernel&#39;s contribution to the Jacobian at the current quadrature point.
void addCoupledVar(const std::string &name, const std::string &doc_string)
This method adds a coupled variable name pair.
virtual void computeJacobian() override
Compute this Kernel&#39;s contribution to the diagonal Jacobian entries.
Definition: Kernel.C:115
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
JacRes
enum to make the code clearer
IntRange< T > make_range(T beg, T end)
static InputParameters validParams()
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
bool isParamValid(const std::string &name) const
Test if the supplied parameter is valid.
Definition: MooseBase.h:199
virtual void computeJacobian() override
Compute this Kernel&#39;s contribution to the diagonal Jacobian entries.