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 30268 : ConservativeAdvectionTempl<is_ad>::generalParams()
21 : {
22 30268 : InputParameters params = GenericKernel<is_ad>::validParams();
23 60536 : 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 121072 : params.addParam<MaterialPropertyName>(
28 : "velocity_scalar_coef",
29 : "1.0",
30 : "Name of material property multiplied against the velocity to scale advection strength.");
31 121072 : MooseEnum upwinding_type("none full", "none");
32 121072 : 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 121072 : params.addParam<MaterialPropertyName>("advected_quantity",
38 : "An optional material property to be advected. If not "
39 : "supplied, then the variable will be used.");
40 90804 : 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 60536 : return params;
45 30268 : }
46 :
47 : template <>
48 : InputParameters
49 15411 : ConservativeAdvectionTempl<false>::validParams()
50 : {
51 15411 : InputParameters params = generalParams();
52 61644 : params.addCoupledVar("velocity", "Velocity vector");
53 92466 : params.deprecateCoupledVar("velocity", "velocity_variable", "12/31/2025");
54 46233 : params.addParam<MaterialPropertyName>("velocity_material", "Velocity vector given as a material");
55 15411 : return params;
56 0 : }
57 :
58 : template <>
59 : InputParameters
60 14857 : ConservativeAdvectionTempl<true>::validParams()
61 : {
62 14857 : InputParameters params = generalParams();
63 59428 : params.addCoupledVar("velocity_variable", "Velocity vector given as a variable");
64 59428 : params.addParam<MaterialPropertyName>("velocity", "Velocity vector given as a material");
65 74285 : params.deprecateParam("velocity", "velocity_material", "12/31/2025");
66 14857 : return params;
67 0 : }
68 :
69 : template <bool is_ad>
70 404 : ConservativeAdvectionTempl<is_ad>::ConservativeAdvectionTempl(const InputParameters & parameters)
71 : : GenericKernel<is_ad>(parameters),
72 404 : _scalar(this->template getGenericMaterialProperty<Real, is_ad>("velocity_scalar_coef")),
73 808 : _coupled_variable_present(isParamValid("velocity_as_variable_gradient")),
74 458 : _coupled_variable_var(_coupled_variable_present ? coupled("velocity_as_variable_gradient") : 0),
75 404 : _velocity(
76 404 : _coupled_variable_present
77 1212 : ? &this->template coupledGenericGradient<is_ad>("velocity_as_variable_gradient")
78 1454 : : (this->isParamValid("velocity_variable")
79 1379 : ? &this->template coupledGenericVectorValue<is_ad>("velocity_variable")
80 629 : : (this->isParamValid("velocity_material")
81 613 : ? &this->template getGenericMaterialProperty<RealVectorValue, is_ad>(
82 : "velocity_material")
83 67 : .get()
84 : : nullptr))),
85 808 : _user_supplied_adv_quant(isParamValid("advected_quantity")),
86 404 : _adv_quant(
87 404 : _user_supplied_adv_quant
88 404 : ? this->template getGenericMaterialProperty<Real, is_ad>("advected_quantity").get()
89 : : _u),
90 404 : _upwinding(
91 808 : this->template getParam<MooseEnum>("upwinding_type").template getEnum<UpwindingType>()),
92 404 : _u_nodal(_var.template genericDofValues<is_ad>()),
93 808 : _upwind_node(0),
94 1616 : _dtotal_mass_out(0)
95 : {
96 404 : if (_coupled_variable_present && _coupled_variable_var == _var.number())
97 8 : 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 684 : 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 892 : if (!_velocity || (_coupled_variable_present && this->isParamValid("velocity_material")) ||
106 1684 : (_coupled_variable_present && this->isParamValid("velocity_variable")) ||
107 2126 : (this->isParamValid("velocity_variable") && this->isParamValid("velocity_material")))
108 16 : 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 392 : 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 392 : }
118 :
119 : template <bool is_ad>
120 : GenericReal<is_ad>
121 19277240 : ConservativeAdvectionTempl<is_ad>::negSpeedQp() const
122 : {
123 19277240 : return -_grad_test[_i][_qp] * (*_velocity)[_qp] * _scalar[_qp];
124 : }
125 :
126 : template <bool is_ad>
127 : GenericReal<is_ad>
128 4753368 : ConservativeAdvectionTempl<is_ad>::computeQpResidual()
129 : {
130 : // This is the no-upwinded version
131 : // It gets called via GenericKernel<is_ad>::computeResidual()
132 4753368 : return negSpeedQp() * _adv_quant[_qp];
133 : }
134 :
135 : template <>
136 : Real
137 9589632 : ConservativeAdvectionTempl<false>::computeQpJacobian()
138 : {
139 : // This is the no-upwinded version
140 : // It gets called via GenericKernel<false>::computeJacobian()
141 9589632 : if (!_user_supplied_adv_quant)
142 9589632 : 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 182944 : ConservativeAdvectionTempl<false>::computeQpOffDiagJacobian(unsigned int jvar)
157 : {
158 : // This is the non-upwinded version
159 : // It gets called via GenericKernel<false>::computeOffDiagJacobian()
160 182944 : if (_coupled_variable_present && _coupled_variable_var == jvar)
161 7840 : return -_grad_test[_i][_qp] * _grad_phi[_j][_qp] * _adv_quant[_qp] * _scalar[_qp];
162 : else
163 175104 : 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 500950 : ConservativeAdvectionTempl<is_ad>::computeResidual()
177 : {
178 500950 : switch (_upwinding)
179 : {
180 339690 : case UpwindingType::none:
181 339690 : GenericKernel<is_ad>::computeResidual();
182 339690 : break;
183 161260 : case UpwindingType::full:
184 161260 : fullUpwind(JacRes::CALCULATE_RESIDUAL);
185 161260 : break;
186 : }
187 500950 : }
188 :
189 : template <bool is_ad>
190 : void
191 314710 : ConservativeAdvectionTempl<is_ad>::computeJacobian()
192 : {
193 314710 : switch (_upwinding)
194 : {
195 154584 : case UpwindingType::none:
196 154584 : GenericKernel<is_ad>::computeJacobian();
197 154584 : break;
198 160126 : case UpwindingType::full:
199 160126 : fullUpwind(JacRes::CALCULATE_JACOBIAN);
200 160126 : break;
201 : }
202 314710 : }
203 :
204 : template <bool is_ad>
205 : void
206 321386 : ConservativeAdvectionTempl<is_ad>::fullUpwind(JacRes res_or_jac)
207 : {
208 : // The number of nodes in the element
209 321386 : 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 321386 : _my_local_re.resize(_var.dofIndices().size());
214 :
215 321386 : if (!is_ad && (res_or_jac == JacRes::CALCULATE_JACOBIAN))
216 160126 : 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 321386 : _upwind_node.resize(num_nodes);
222 1567258 : for (_i = 0; _i < num_nodes; ++_i)
223 : {
224 6180112 : for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
225 4934240 : _my_local_re(_i) += this->_JxW[_qp] * this->_coord[_qp] * negSpeedQp();
226 1245872 : _upwind_node[_i] = (MetaPhysicL::raw_value(_my_local_re(_i)) >= 0.0);
227 : }
228 :
229 : // Variables used to ensure mass conservation
230 321386 : GenericReal<is_ad> total_mass_out = 0.0;
231 321386 : GenericReal<is_ad> total_in = 0.0;
232 321386 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
233 160126 : _dtotal_mass_out.assign(num_nodes, 0.0);
234 :
235 1567258 : for (const auto n : make_range(num_nodes))
236 : {
237 1245872 : if (_upwind_node[n])
238 : {
239 : if constexpr (!is_ad)
240 623080 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
241 : {
242 309272 : 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 309272 : _local_ke(n, n) += _my_local_re(n);
250 :
251 309272 : _dtotal_mass_out[n] += _local_ke(n, n);
252 : }
253 623080 : _my_local_re(n) *= _u_nodal[n];
254 623080 : total_mass_out += _my_local_re(n);
255 : }
256 : else // downwind node
257 622792 : 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 1567258 : for (const auto n : make_range(num_nodes))
263 1245872 : if (!_upwind_node[n]) // downwind node
264 : {
265 : if constexpr (!is_ad)
266 622792 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
267 1524472 : for (_j = 0; _j < _phi.size(); _j++)
268 1215344 : _local_ke(n, _j) += _my_local_re(n) * _dtotal_mass_out[_j] / total_in;
269 622792 : _my_local_re(n) *= total_mass_out / total_in;
270 : }
271 :
272 : // Add the result to the residual and jacobian
273 321386 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
274 : {
275 161260 : this->addResiduals(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
276 :
277 161260 : if (this->_has_save_in)
278 : {
279 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
280 0 : for (const auto & var : this->_save_in)
281 0 : var->sys().solution().add_vector(MetaPhysicL::raw_value(_my_local_re), var->dofIndices());
282 0 : }
283 : }
284 :
285 321386 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
286 : {
287 : if constexpr (is_ad)
288 0 : this->addJacobian(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
289 : else
290 160126 : accumulateTaggedLocalMatrix();
291 : }
292 321386 : }
293 :
294 : template class ConservativeAdvectionTempl<false>;
295 : template class ConservativeAdvectionTempl<true>;
|
