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 "ConservativeAdvection.h"
11 : #include "SystemBase.h"
12 :
13 : using namespace libMesh;
14 :
15 : registerMooseObject("MooseApp", ConservativeAdvection);
16 : registerMooseObject("MooseApp", ADConservativeAdvection);
17 :
18 : template <bool is_ad>
19 : InputParameters
20 6856 : ConservativeAdvectionTempl<is_ad>::generalParams()
21 : {
22 6856 : InputParameters params = GenericKernel<is_ad>::validParams();
23 13712 : 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 27424 : params.addParam<MaterialPropertyName>(
28 : "velocity_scalar_coef",
29 : "1.0",
30 : "Name of material property multiplied against the velocity to scale advection strength.");
31 27424 : MooseEnum upwinding_type("none full", "none");
32 27424 : 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 27424 : params.addParam<MaterialPropertyName>("advected_quantity",
38 : "An optional material property to be advected. If not "
39 : "supplied, then the variable will be used.");
40 20568 : 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 13712 : return params;
45 6856 : }
46 :
47 : template <>
48 : InputParameters
49 3695 : ConservativeAdvectionTempl<false>::validParams()
50 : {
51 3695 : InputParameters params = generalParams();
52 14780 : params.addCoupledVar("velocity", "Velocity vector");
53 22170 : params.deprecateCoupledVar("velocity", "velocity_variable", "12/31/2025");
54 11085 : params.addParam<MaterialPropertyName>("velocity_material", "Velocity vector given as a material");
55 3695 : return params;
56 0 : }
57 :
58 : template <>
59 : InputParameters
60 3161 : ConservativeAdvectionTempl<true>::validParams()
61 : {
62 3161 : InputParameters params = generalParams();
63 12644 : params.addCoupledVar("velocity_variable", "Velocity vector given as a variable");
64 12644 : params.addParam<MaterialPropertyName>("velocity", "Velocity vector given as a material");
65 15805 : params.deprecateParam("velocity", "velocity_material", "12/31/2025");
66 3161 : return params;
67 0 : }
68 :
69 : template <bool is_ad>
70 384 : ConservativeAdvectionTempl<is_ad>::ConservativeAdvectionTempl(const InputParameters & parameters)
71 : : GenericKernel<is_ad>(parameters),
72 384 : _scalar(this->template getGenericMaterialProperty<Real, is_ad>("velocity_scalar_coef")),
73 768 : _coupled_variable_present(isParamValid("velocity_as_variable_gradient")),
74 433 : _coupled_variable_var(_coupled_variable_present ? coupled("velocity_as_variable_gradient") : 0),
75 384 : _velocity(
76 384 : _coupled_variable_present
77 1152 : ? &this->template coupledGenericGradient<is_ad>("velocity_as_variable_gradient")
78 1389 : : (this->isParamValid("velocity_variable")
79 1321 : ? &this->template coupledGenericVectorValue<is_ad>("velocity_variable")
80 588 : : (this->isParamValid("velocity_material")
81 576 : ? &this->template getGenericMaterialProperty<RealVectorValue, is_ad>(
82 : "velocity_material")
83 62 : .get()
84 : : nullptr))),
85 768 : _user_supplied_adv_quant(isParamValid("advected_quantity")),
86 384 : _adv_quant(
87 384 : _user_supplied_adv_quant
88 384 : ? this->template getGenericMaterialProperty<Real, is_ad>("advected_quantity").get()
89 : : _u),
90 384 : _upwinding(
91 768 : this->template getParam<MooseEnum>("upwinding_type").template getEnum<UpwindingType>()),
92 384 : _u_nodal(_var.template genericDofValues<is_ad>()),
93 768 : _upwind_node(0),
94 1536 : _dtotal_mass_out(0)
95 : {
96 384 : if (_coupled_variable_present && _coupled_variable_var == _var.number())
97 6 : 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 643 : if (_upwinding != UpwindingType::none && this->isParamValid("advected_quantity"))
101 0 : paramError(
102 : "advected_quantity",
103 : "Upwinding is not compatible with an advected quantity that is not the primary variable.");
104 :
105 848 : if (!_velocity || (_coupled_variable_present && this->isParamValid("velocity_material")) ||
106 1604 : (_coupled_variable_present && this->isParamValid("velocity_variable")) ||
107 2040 : (this->isParamValid("velocity_variable") && this->isParamValid("velocity_material")))
108 12 : 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 375 : if (this->_has_diag_save_in)
114 0 : 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 375 : }
118 :
119 : template <bool is_ad>
120 : GenericReal<is_ad>
121 29543744 : ConservativeAdvectionTempl<is_ad>::negSpeedQp() const
122 : {
123 29543744 : return -_grad_test[_i][_qp] * (*_velocity)[_qp] * _scalar[_qp];
124 : }
125 :
126 : template <bool is_ad>
127 : GenericReal<is_ad>
128 8324648 : ConservativeAdvectionTempl<is_ad>::computeQpResidual()
129 : {
130 : // This is the no-upwinded version
131 : // It gets called via GenericKernel<is_ad>::computeResidual()
132 8324648 : return negSpeedQp() * _adv_quant[_qp];
133 : }
134 :
135 : template <>
136 : Real
137 16829184 : ConservativeAdvectionTempl<false>::computeQpJacobian()
138 : {
139 : // This is the no-upwinded version
140 : // It gets called via GenericKernel<false>::computeJacobian()
141 16829184 : if (!_user_supplied_adv_quant)
142 16829184 : return negSpeedQp() * _phi[_j][_qp];
143 0 : return 0.0;
144 : }
145 :
146 : template <>
147 : Real
148 0 : ConservativeAdvectionTempl<true>::computeQpJacobian()
149 : {
150 0 : mooseError("Internal error, should never get here when using AD");
151 : return 0.0;
152 : }
153 :
154 : template <>
155 : Real
156 163632 : ConservativeAdvectionTempl<false>::computeQpOffDiagJacobian(unsigned int jvar)
157 : {
158 : // This is the non-upwinded version
159 : // It gets called via GenericKernel<false>::computeOffDiagJacobian()
160 163632 : if (_coupled_variable_present && _coupled_variable_var == jvar)
161 6960 : return -_grad_test[_i][_qp] * _grad_phi[_j][_qp] * _adv_quant[_qp] * _scalar[_qp];
162 : else
163 156672 : return 0.0;
164 : }
165 :
166 : template <>
167 : Real
168 0 : ConservativeAdvectionTempl<true>::computeQpOffDiagJacobian(unsigned int /*jvar*/)
169 : {
170 0 : mooseError("Internal error, should never get here when using AD");
171 : return 0.0;
172 : }
173 :
174 : template <bool is_ad>
175 : void
176 700570 : ConservativeAdvectionTempl<is_ad>::computeResidual()
177 : {
178 700570 : switch (_upwinding)
179 : {
180 557074 : case UpwindingType::none:
181 557074 : GenericKernel<is_ad>::computeResidual();
182 557074 : break;
183 143496 : case UpwindingType::full:
184 143496 : fullUpwind(JacRes::CALCULATE_RESIDUAL);
185 143496 : break;
186 : }
187 700570 : }
188 :
189 : template <bool is_ad>
190 : void
191 409708 : ConservativeAdvectionTempl<is_ad>::computeJacobian()
192 : {
193 409708 : switch (_upwinding)
194 : {
195 267184 : case UpwindingType::none:
196 267184 : GenericKernel<is_ad>::computeJacobian();
197 267184 : break;
198 142524 : case UpwindingType::full:
199 142524 : fullUpwind(JacRes::CALCULATE_JACOBIAN);
200 142524 : break;
201 : }
202 409708 : }
203 :
204 : template <bool is_ad>
205 : void
206 286020 : ConservativeAdvectionTempl<is_ad>::fullUpwind(JacRes res_or_jac)
207 : {
208 : // The number of nodes in the element
209 286020 : 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 286020 : _my_local_re.resize(_var.dofIndices().size());
214 :
215 286020 : if (!is_ad && (res_or_jac == JacRes::CALCULATE_JACOBIAN))
216 142524 : 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 286020 : _upwind_node.resize(num_nodes);
222 1394748 : for (_i = 0; _i < num_nodes; ++_i)
223 : {
224 5498640 : for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
225 4389912 : _my_local_re(_i) += this->_JxW[_qp] * this->_coord[_qp] * negSpeedQp();
226 1108728 : _upwind_node[_i] = (MetaPhysicL::raw_value(_my_local_re(_i)) >= 0.0);
227 : }
228 :
229 : // Variables used to ensure mass conservation
230 286020 : GenericReal<is_ad> total_mass_out = 0.0;
231 286020 : GenericReal<is_ad> total_in = 0.0;
232 286020 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
233 142524 : _dtotal_mass_out.assign(num_nodes, 0.0);
234 :
235 1394748 : for (const auto n : make_range(num_nodes))
236 : {
237 1108728 : if (_upwind_node[n])
238 : {
239 : if constexpr (!is_ad)
240 554472 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
241 : {
242 275292 : 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 275292 : _local_ke(n, n) += _my_local_re(n);
250 :
251 275292 : _dtotal_mass_out[n] += _local_ke(n, n);
252 : }
253 554472 : _my_local_re(n) *= getUNodal(n);
254 554472 : total_mass_out += _my_local_re(n);
255 : }
256 : else // downwind node
257 554256 : 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 1394748 : for (const auto n : make_range(num_nodes))
263 1108728 : if (!_upwind_node[n]) // downwind node
264 : {
265 : if constexpr (!is_ad)
266 554256 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
267 1357056 : for (_j = 0; _j < _phi.size(); _j++)
268 1081872 : _local_ke(n, _j) += _my_local_re(n) * _dtotal_mass_out[_j] / total_in;
269 554256 : _my_local_re(n) *= total_mass_out / total_in;
270 : }
271 :
272 : // Add the result to the residual and jacobian
273 286020 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
274 : {
275 143496 : this->addResiduals(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
276 :
277 143496 : if (this->_has_save_in)
278 0 : for (const auto & var : this->_save_in)
279 0 : var->sys().solution().add_vector(MetaPhysicL::raw_value(_my_local_re), var->dofIndices());
280 : }
281 :
282 286020 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
283 : {
284 : if constexpr (is_ad)
285 0 : this->addJacobian(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
286 : else
287 142524 : accumulateTaggedLocalMatrix();
288 : }
289 286020 : }
290 :
291 : template class ConservativeAdvectionTempl<false>;
292 : template class ConservativeAdvectionTempl<true>;
|
