www.mooseframework.org
ConservativeAdvection.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
13 
14 template <>
17 {
19  params.addClassDescription("Conservative form of $\\nabla \\cdot \\vec{v} u$ which in its weak "
20  "form is given by: $(-\\nabla \\psi_i, \\vec{v} u)$.");
21  params.addRequiredParam<RealVectorValue>("velocity", "Velocity vector");
22  MooseEnum upwinding_type("none full", "none");
23  params.addParam<MooseEnum>("upwinding_type",
24  upwinding_type,
25  "Type of upwinding used. None: Typically results in overshoots and "
26  "undershoots, but numerical diffusion is minimized. Full: Overshoots "
27  "and undershoots are avoided, but numerical diffusion is large");
28  return params;
29 }
30 
32  : Kernel(parameters),
33  _velocity(getParam<RealVectorValue>("velocity")),
34  _upwinding(getParam<MooseEnum>("upwinding_type").getEnum<UpwindingType>()),
35  _u_nodal(_var.dofValues()),
36  _upwind_node(0),
37  _dtotal_mass_out(0)
38 {
39 }
40 
41 Real
43 {
44  return -_grad_test[_i][_qp] * _velocity;
45 }
46 
47 Real
49 {
50  // This is the no-upwinded version
51  // It gets called via Kernel::computeResidual()
52  return negSpeedQp() * _u[_qp];
53 }
54 
55 Real
57 {
58  // This is the no-upwinded version
59  // It gets called via Kernel::computeJacobian()
60  return negSpeedQp() * _phi[_j][_qp];
61 }
62 
63 void
65 {
66  switch (_upwinding)
67  {
70  break;
73  break;
74  }
75 }
76 
77 void
79 {
80  switch (_upwinding)
81  {
84  break;
87  break;
88  }
89 }
90 
91 void
93 {
94  // The number of nodes in the element
95  const unsigned int num_nodes = _test.size();
96 
97  // Even if we are computing the Jacobian we still need to compute the outflow from each node to
98  // see which nodes are upwind and which are downwind
100 
101  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
103 
104  // Compute the outflux from each node and store in _local_re
105  // If _local_re is positive at the node, mass (or whatever the Variable represents) is flowing out
106  // of the node
107  _upwind_node.resize(num_nodes);
108  for (_i = 0; _i < num_nodes; ++_i)
109  {
110  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
111  _local_re(_i) += _JxW[_qp] * _coord[_qp] * negSpeedQp();
112  _upwind_node[_i] = (_local_re(_i) >= 0.0);
113  }
114 
115  // Variables used to ensure mass conservation
116  Real total_mass_out = 0.0;
117  Real total_in = 0.0;
118  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
119  _dtotal_mass_out.assign(num_nodes, 0.0);
120 
121  for (unsigned int n = 0; n < num_nodes; ++n)
122  {
123  if (_upwind_node[n])
124  {
125  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
126  {
127  if (_test.size() == _phi.size())
128  /* u at node=n depends only on the u at node=n, by construction. For
129  * linear-lagrange variables, this means that Jacobian entries involving the derivative
130  * will only be nonzero for derivatives wrt variable at node=n. Hence the
131  * (n, n) in the line below. The above "if" statement catches other variable types
132  * (eg constant monomials)
133  */
134  _local_ke(n, n) += _local_re(n);
135 
137  }
138  _local_re(n) *= _u_nodal[n];
139  total_mass_out += _local_re(n);
140  }
141  else // downwind node
142  total_in -= _local_re(n); // note the -= means the result is positive
143  }
144 
145  // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
146  // weighted by their local_re values
147  for (unsigned int n = 0; n < num_nodes; ++n)
148  {
149  if (!_upwind_node[n]) // downwind node
150  {
151  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
152  for (_j = 0; _j < _phi.size(); _j++)
153  _local_ke(n, _j) += _local_re(n) * _dtotal_mass_out[_j] / total_in;
154  _local_re(n) *= total_mass_out / total_in;
155  }
156  }
157 
158  // Add the result to the residual and jacobian
159  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
160  {
162 
163  if (_has_save_in)
164  {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  for (const auto & var : _save_in)
167  var->sys().solution().add_vector(_local_re, var->dofIndices());
168  }
169  }
170 
171  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
172  {
174 
175  if (_has_diag_save_in)
176  {
177  unsigned int rows = _local_ke.m();
178  DenseVector<Number> diag(rows);
179  for (unsigned int i = 0; i < rows; i++)
180  diag(i) = _local_ke(i, i);
181 
182  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
183  for (const auto & var : _diag_save_in)
184  var->sys().solution().add_vector(diag, var->dofIndices());
185  }
186  }
187 }
RealVectorValue _velocity
advection velocity
virtual void computeResidual() override
Compute this Kernel&#39;s contribution to the residual.
Definition: Kernel.C:92
virtual void computeResidual() override
Compute this Kernel&#39;s contribution to the residual.
void accumulateTaggedLocalResidual()
Local residual blocks will be appended by adding the current local kernel residual.
MooseVariable & _var
This is a regular kernel so we cast to a regular MooseVariable.
Definition: Kernel.h:54
std::vector< MooseVariableFEBase * > _save_in
Definition: KernelBase.h:172
const MooseArray< Real > & _JxW
The current quadrature point weight value.
Definition: KernelBase.h:159
VectorValue< Real > RealVectorValue
Definition: Assembly.h:31
unsigned int number() const
Get variable number coming from libMesh.
enum ConservativeAdvection::UpwindingType _upwinding
virtual Real computeQpResidual() override
Compute this Kernel&#39;s contribution to the residual at the current quadrature point.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
ConservativeAdvection(const InputParameters &parameters)
bool _has_diag_save_in
The aux variables to save the diagonal Jacobian contributions to.
Definition: KernelBase.h:176
DenseMatrix< Number > _local_ke
Holds residual entries as they are accumulated by this Kernel.
virtual Real computeQpJacobian() override
Compute this Kernel&#39;s contribution to the Jacobian at the current quadrature point.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
registerMooseObject("MooseApp", ConservativeAdvection)
Assembly & _assembly
Reference to this Kernel&#39;s assembly object.
Definition: KernelBase.h:139
const VariableTestValue & _test
the current test function
Definition: Kernel.h:57
virtual void computeJacobian() override
Compute this Kernel&#39;s contribution to the diagonal Jacobian entries.
std::vector< MooseVariableFEBase * > _diag_save_in
Definition: KernelBase.h:177
void fullUpwind(JacRes res_or_jac)
Calculates the fully-upwind Residual and Jacobian (depending on res_or_jac)
std::vector< Real > _dtotal_mass_out
In the full-upwind scheme d(total_mass_out)/d(variable_at_node_i)
const QBase *const & _qrule
active quadrature rule
Definition: KernelBase.h:156
bool _has_save_in
The aux variables to save the residual contributions to.
Definition: KernelBase.h:171
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
Real negSpeedQp() const
Returns - _grad_test * velocity.
unsigned int _i
current index for the test function
Definition: KernelBase.h:165
void accumulateTaggedLocalMatrix()
Local Jacobian blocks will be appended by adding the current local kernel Jacobian.
const VariableValue & _u_nodal
Nodal value of u, used for full upwinding.
const MooseArray< Real > & _coord
The scaling factor to convert from cartesian to another coordinate system (e.g rz, spherical, etc.)
Definition: KernelBase.h:162
UpwindingType
Type of upwinding.
virtual void computeJacobian() override
Compute this Kernel&#39;s contribution to the diagonal Jacobian entries.
Definition: Kernel.C:112
unsigned int _j
current index for the shape function
Definition: KernelBase.h:168
InputParameters validParams< Kernel >()
Definition: Kernel.C:24
std::vector< bool > _upwind_node
In the full-upwind scheme, whether a node is an upwind node.
const VariableTestGradient & _grad_test
gradient of the test function
Definition: Kernel.h:60
PetscInt n
DenseVector< Number > _local_re
Holds residual entries as they are accumulated by this Kernel.
Definition: Kernel.h:20
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 option parameter and a documentation string to the InputParameters object...
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
Prepare data for computing element residual the according to active tags.
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
Prepare data for computing element jacobian according to the ative tags.
const VariablePhiValue & _phi
the current shape functions
Definition: Kernel.h:63
InputParameters validParams< ConservativeAdvection >()
Advection of the variable by the velocity provided by the user.
const VariableValue & _u
Holds the solution at current quadrature points.
Definition: Kernel.h:69
unsigned int _qp
The current quadrature point index.
Definition: KernelBase.h:150
JacRes
enum to make the code clearer