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 29115 : ConservativeAdvectionTempl<is_ad>::generalParams()
21 : {
22 29115 : InputParameters params = GenericKernel<is_ad>::validParams();
23 29115 : 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)$.");
25 29115 : MooseEnum upwinding_type("none full", "none");
26 29115 : params.addParam<MooseEnum>("upwinding_type",
27 : upwinding_type,
28 : "Type of upwinding used. None: Typically results in overshoots and "
29 : "undershoots, but numerical diffusion is minimized. Full: Overshoots "
30 : "and undershoots are avoided, but numerical diffusion is large");
31 29115 : params.addParam<MaterialPropertyName>("advected_quantity",
32 : "An optional material property to be advected. If not "
33 : "supplied, then the variable will be used.");
34 58230 : return params;
35 29115 : }
36 :
37 : template <>
38 : InputParameters
39 14767 : ConservativeAdvectionTempl<false>::validParams()
40 : {
41 14767 : InputParameters params = generalParams();
42 14767 : params.addCoupledVar("velocity", "Velocity vector");
43 14767 : params.deprecateCoupledVar("velocity", "velocity_variable", "12/31/2025");
44 14767 : params.addParam<MaterialPropertyName>("velocity_material", "Velocity vector given as a material");
45 14767 : return params;
46 0 : }
47 :
48 : template <>
49 : InputParameters
50 14348 : ConservativeAdvectionTempl<true>::validParams()
51 : {
52 14348 : InputParameters params = generalParams();
53 14348 : params.addCoupledVar("velocity_variable", "Velocity vector given as a variable");
54 14348 : params.addParam<MaterialPropertyName>("velocity", "Velocity vector given as a material");
55 14348 : params.deprecateParam("velocity", "velocity_material", "12/31/2025");
56 14348 : return params;
57 0 : }
58 :
59 : template <bool is_ad>
60 304 : ConservativeAdvectionTempl<is_ad>::ConservativeAdvectionTempl(const InputParameters & parameters)
61 : : GenericKernel<is_ad>(parameters),
62 608 : _velocity(this->isParamValid("velocity_variable")
63 655 : ? &this->template coupledGenericVectorValue<is_ad>("velocity_variable")
64 351 : : (this->isParamValid("velocity_material")
65 351 : ? &this->template getGenericMaterialProperty<RealVectorValue, is_ad>(
66 : "velocity_material")
67 39 : .get()
68 : : nullptr)),
69 304 : _adv_quant(
70 608 : isParamValid("advected_quantity")
71 304 : ? this->template getGenericMaterialProperty<Real, is_ad>("advected_quantity").get()
72 : : _u),
73 304 : _upwinding(
74 304 : this->template getParam<MooseEnum>("upwinding_type").template getEnum<UpwindingType>()),
75 304 : _u_nodal(_var.template genericDofValues<is_ad>()),
76 304 : _upwind_node(0),
77 608 : _dtotal_mass_out(0)
78 : {
79 304 : if (_upwinding != UpwindingType::none && this->isParamValid("advected_quantity"))
80 0 : paramError(
81 : "advected_quantity",
82 : "Upwinding is not compatable with an advected quantity that is not the primary variable.");
83 :
84 600 : if (!_velocity ||
85 600 : (this->isParamValid("velocity_variable") && this->isParamValid("velocity_material")))
86 8 : paramError("velocity_variable",
87 : "Either velocity_variable or velocity_material must be specified");
88 :
89 296 : if (this->_has_diag_save_in)
90 0 : paramError("diag_save_in",
91 : "_local_ke not computed for global AD indexing. Save-in is deprecated anyway. Use "
92 : "the tagging system instead.");
93 296 : }
94 :
95 : template <bool is_ad>
96 : GenericReal<is_ad>
97 15400216 : ConservativeAdvectionTempl<is_ad>::negSpeedQp() const
98 : {
99 15400216 : return -_grad_test[_i][_qp] * (*_velocity)[_qp];
100 : }
101 :
102 : template <bool is_ad>
103 : GenericReal<is_ad>
104 2704704 : ConservativeAdvectionTempl<is_ad>::computeQpResidual()
105 : {
106 : // This is the no-upwinded version
107 : // It gets called via GenericKernel<is_ad>::computeResidual()
108 2704704 : return negSpeedQp() * _adv_quant[_qp];
109 : }
110 :
111 : template <>
112 : Real
113 8352384 : ConservativeAdvectionTempl<false>::computeQpJacobian()
114 : {
115 : // This is the no-upwinded version
116 : // It gets called via GenericKernel<false>::computeJacobian()
117 8352384 : return negSpeedQp() * _phi[_j][_qp];
118 : }
119 :
120 : template <>
121 : Real
122 0 : ConservativeAdvectionTempl<true>::computeQpJacobian()
123 : {
124 0 : mooseError("Internal error, should never get here when using AD");
125 : return 0.0;
126 : }
127 :
128 : template <bool is_ad>
129 : void
130 314480 : ConservativeAdvectionTempl<is_ad>::computeResidual()
131 : {
132 314480 : switch (_upwinding)
133 : {
134 172716 : case UpwindingType::none:
135 172716 : GenericKernel<is_ad>::computeResidual();
136 172716 : break;
137 141764 : case UpwindingType::full:
138 141764 : fullUpwind(JacRes::CALCULATE_RESIDUAL);
139 141764 : break;
140 : }
141 314480 : }
142 :
143 : template <bool is_ad>
144 : void
145 272942 : ConservativeAdvectionTempl<is_ad>::computeJacobian()
146 : {
147 272942 : switch (_upwinding)
148 : {
149 132312 : case UpwindingType::none:
150 132312 : GenericKernel<is_ad>::computeJacobian();
151 132312 : break;
152 140630 : case UpwindingType::full:
153 140630 : fullUpwind(JacRes::CALCULATE_JACOBIAN);
154 140630 : break;
155 : }
156 272942 : }
157 :
158 : template <bool is_ad>
159 : void
160 282394 : ConservativeAdvectionTempl<is_ad>::fullUpwind(JacRes res_or_jac)
161 : {
162 : // The number of nodes in the element
163 282394 : const unsigned int num_nodes = _test.size();
164 :
165 : // Even if we are computing the Jacobian we still need to compute the outflow from each node to
166 : // see which nodes are upwind and which are downwind
167 282394 : _my_local_re.resize(_var.dofIndices().size());
168 :
169 282394 : if (!is_ad && (res_or_jac == JacRes::CALCULATE_JACOBIAN))
170 140630 : prepareMatrixTag(this->_assembly, _var.number(), _var.number());
171 :
172 : // Compute the outflux from each node and store in _my_local_re
173 : // If _my_local_re is positive at the node, mass (or whatever the Variable represents) is flowing
174 : // out of the node
175 282394 : _upwind_node.resize(num_nodes);
176 1377794 : for (_i = 0; _i < num_nodes; ++_i)
177 : {
178 5438528 : for (_qp = 0; _qp < this->_qrule->n_points(); _qp++)
179 4343128 : _my_local_re(_i) += this->_JxW[_qp] * this->_coord[_qp] * negSpeedQp();
180 1095400 : _upwind_node[_i] = (MetaPhysicL::raw_value(_my_local_re(_i)) >= 0.0);
181 : }
182 :
183 : // Variables used to ensure mass conservation
184 282394 : GenericReal<is_ad> total_mass_out = 0.0;
185 282394 : GenericReal<is_ad> total_in = 0.0;
186 282394 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
187 140630 : _dtotal_mass_out.assign(num_nodes, 0.0);
188 :
189 1377794 : for (const auto n : make_range(num_nodes))
190 : {
191 1095400 : if (_upwind_node[n])
192 : {
193 : if constexpr (!is_ad)
194 547808 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
195 : {
196 271636 : if (_test.size() == _phi.size())
197 : /* u at node=n depends only on the u at node=n, by construction. For
198 : * linear-lagrange variables, this means that Jacobian entries involving the derivative
199 : * will only be nonzero for derivatives wrt variable at node=n. Hence the
200 : * (n, n) in the line below. The above "if" statement catches other variable types
201 : * (eg constant monomials)
202 : */
203 271636 : _local_ke(n, n) += _my_local_re(n);
204 :
205 271636 : _dtotal_mass_out[n] += _local_ke(n, n);
206 : }
207 547808 : _my_local_re(n) *= _u_nodal[n];
208 547808 : total_mass_out += _my_local_re(n);
209 : }
210 : else // downwind node
211 547592 : total_in -= _my_local_re(n); // note the -= means the result is positive
212 : }
213 :
214 : // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
215 : // weighted by their local_re values
216 1377794 : for (const auto n : make_range(num_nodes))
217 1095400 : if (!_upwind_node[n]) // downwind node
218 : {
219 : if constexpr (!is_ad)
220 547592 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
221 1339112 : for (_j = 0; _j < _phi.size(); _j++)
222 1067584 : _local_ke(n, _j) += _my_local_re(n) * _dtotal_mass_out[_j] / total_in;
223 547592 : _my_local_re(n) *= total_mass_out / total_in;
224 : }
225 :
226 : // Add the result to the residual and jacobian
227 282394 : if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
228 : {
229 141764 : this->addResiduals(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
230 :
231 141764 : if (this->_has_save_in)
232 : {
233 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
234 0 : for (const auto & var : this->_save_in)
235 0 : var->sys().solution().add_vector(MetaPhysicL::raw_value(_my_local_re), var->dofIndices());
236 0 : }
237 : }
238 :
239 282394 : if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
240 : {
241 : if constexpr (is_ad)
242 0 : this->addJacobian(this->_assembly, _my_local_re, _var.dofIndices(), _var.scalingFactor());
243 : else
244 140630 : accumulateTaggedLocalMatrix();
245 : }
246 282394 : }
247 :
248 : template class ConservativeAdvectionTempl<false>;
249 : template class ConservativeAdvectionTempl<true>;
|
