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 "INSFVMomentumAdvection.h"
11 : #include "NS.h"
12 : #include "FVUtils.h"
13 : #include "INSFVRhieChowInterpolator.h"
14 : #include "SystemBase.h"
15 : #include "NSFVUtils.h"
16 :
17 : registerMooseObject("NavierStokesApp", INSFVMomentumAdvection);
18 :
19 : InputParameters
20 79316 : INSFVMomentumAdvection::uniqueParams()
21 : {
22 79316 : auto params = emptyInputParameters();
23 158632 : params.addParam<Real>(
24 : "characteristic_speed",
25 : "The characteristic speed. For porous medium simulations, this characteristic speed should "
26 : "correspond to the superficial velocity, not the interstitial.");
27 79316 : return params;
28 0 : }
29 :
30 : std::vector<std::string>
31 1621 : INSFVMomentumAdvection::listOfCommonParams()
32 : {
33 1621 : return {"characteristic_speed"};
34 : }
35 :
36 : InputParameters
37 32294 : INSFVMomentumAdvection::validParams()
38 : {
39 32294 : InputParameters params = INSFVAdvectionKernel::validParams();
40 32294 : params += INSFVMomentumResidualObject::validParams();
41 32294 : params += INSFVMomentumAdvection::uniqueParams();
42 32294 : params.addRequiredParam<MooseFunctorName>(NS::density, "Density functor");
43 32294 : params.addClassDescription("Object for advecting momentum, e.g. rho*u");
44 32294 : return params;
45 0 : }
46 :
47 16938 : INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
48 : : INSFVAdvectionKernel(params),
49 : INSFVMomentumResidualObject(*this),
50 33876 : _rho(getFunctor<ADReal>(NS::density)),
51 33876 : _approximate_as(isParamValid("characteristic_speed")),
52 35084 : _cs(_approximate_as ? getParam<Real>("characteristic_speed") : 0)
53 : {
54 16938 : }
55 :
56 : void
57 22644 : INSFVMomentumAdvection::initialSetup()
58 : {
59 22644 : INSFVAdvectionKernel::initialSetup();
60 :
61 22644 : _rho_u = std::make_unique<PiecewiseByBlockLambdaFunctor<ADReal>>(
62 : "rho_u",
63 386123466 : [this](const auto & r, const auto & t) -> ADReal
64 386078178 : { return _rho(r, t) * _var(r, t) / epsilon()(r, t); },
65 67932 : std::set<ExecFlagType>({EXEC_ALWAYS}),
66 : _mesh,
67 22644 : this->blockIDs());
68 45288 : }
69 :
70 : ADReal
71 0 : INSFVMomentumAdvection::computeQpResidual()
72 : {
73 0 : mooseError("Should never be called");
74 : }
75 :
76 : void
77 190408499 : INSFVMomentumAdvection::computeResidualsAndAData(const FaceInfo & fi)
78 : {
79 : mooseAssert(!skipForBoundary(fi),
80 : "We shouldn't get in here if we're supposed to skip for a boundary");
81 :
82 190408499 : _face_info = &fi;
83 190408499 : _normal = fi.normal();
84 190408499 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
85 190408499 : const auto state = determineState();
86 :
87 : using namespace Moose::FV;
88 :
89 190408499 : const bool correct_skewness = _advected_interp_method == InterpMethod::SkewCorrectedAverage;
90 190408499 : const bool on_boundary = onBoundary(fi);
91 :
92 190408499 : _elem_residual = 0, _neighbor_residual = 0, _ae = 0, _an = 0;
93 :
94 190408499 : const auto v_face = velocity();
95 190408499 : const auto vdotn = _normal * v_face;
96 :
97 190408499 : if (on_boundary)
98 : {
99 3249594 : const auto ssf = singleSidedFaceArg();
100 3249594 : const Elem * const sided_elem = ssf.face_side;
101 3249594 : const auto dof_number = sided_elem->dof_number(_sys.number(), _var.number(), 0);
102 3249594 : const auto rho_face = _rho(ssf, state);
103 3249594 : const auto eps_face = epsilon()(ssf, state);
104 3249594 : const auto u_face = _var(ssf, state);
105 3249594 : const Real d_u_face_d_dof = u_face.derivatives()[dof_number];
106 3249594 : const auto coeff = vdotn * rho_face / eps_face;
107 :
108 3249594 : if (sided_elem == fi.elemPtr())
109 : {
110 3249094 : _ae = coeff * d_u_face_d_dof;
111 3249094 : _elem_residual = coeff * u_face;
112 3249094 : if (_approximate_as)
113 12960 : _ae = _cs / fi.elem().n_sides() * rho_face / eps_face;
114 : }
115 : else
116 : {
117 1000 : _an = -coeff * d_u_face_d_dof;
118 1000 : _neighbor_residual = -coeff * u_face;
119 500 : if (_approximate_as)
120 0 : _an = _cs / fi.neighborPtr()->n_sides() * rho_face / eps_face;
121 : }
122 : }
123 : else // we are an internal fluid flow face
124 : {
125 187158905 : const bool elem_is_upwind = MetaPhysicL::raw_value(v_face) * _normal > 0;
126 187158905 : const auto & limiter_time = _subproblem.isTransient()
127 187158905 : ? Moose::StateArg(1, Moose::SolutionIterationType::Time)
128 : : Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear);
129 187158905 : const Moose::FaceArg advected_face_arg{&fi,
130 187158905 : limiterType(_advected_interp_method),
131 : elem_is_upwind,
132 : correct_skewness,
133 : nullptr,
134 187158905 : &limiter_time};
135 187158905 : if (const auto [is_jump, eps_elem_face, eps_neighbor_face] =
136 187158905 : NS::isPorosityJumpFace(epsilon(), fi, state);
137 187158905 : is_jump)
138 : {
139 : // For a weakly compressible formulation, the density should not depend on pressure and
140 : // consequently the density should not be impacted by the pressure jump that occurs at a
141 : // porosity jump. Consequently we will allow evaluation of the density using both upstream and
142 : // downstream information
143 31176 : const auto rho_face = _rho(advected_face_arg, state);
144 :
145 : // We set the + and - sides of the superficial velocity equal to the interpolated value
146 31176 : const auto & var_elem_face = v_face(_index);
147 : const auto & var_neighbor_face = v_face(_index);
148 :
149 31176 : const auto elem_dof_number = fi.elem().dof_number(_sys.number(), _var.number(), 0);
150 31176 : const auto neighbor_dof_number = fi.neighbor().dof_number(_sys.number(), _var.number(), 0);
151 :
152 62352 : const auto d_var_elem_face_d_elem_dof = var_elem_face.derivatives()[elem_dof_number];
153 : const auto d_var_neighbor_face_d_neighbor_dof =
154 31176 : var_neighbor_face.derivatives()[neighbor_dof_number];
155 :
156 : // We allow a discontintuity in the advective momentum flux at the jump face. Mainly the
157 : // advective flux is:
158 : // elem side:
159 : // rho * v_superficial / eps_elem * v_superficial = rho * v_interstitial_elem * v_superficial
160 : // neighbor side:
161 : // rho * v_superficial / eps_neigh * v_superficial = rho * v_interstitial_neigh *
162 : // v_superficial
163 31176 : const auto elem_coeff = vdotn * rho_face / eps_elem_face;
164 31176 : const auto neighbor_coeff = vdotn * rho_face / eps_neighbor_face;
165 31176 : _ae = elem_coeff * d_var_elem_face_d_elem_dof;
166 31176 : _elem_residual = elem_coeff * var_elem_face;
167 62352 : _an = -neighbor_coeff * d_var_neighbor_face_d_neighbor_dof;
168 62352 : _neighbor_residual = -neighbor_coeff * var_neighbor_face;
169 31176 : if (_approximate_as)
170 : {
171 0 : _ae = _cs / fi.elem().n_sides() * rho_face / eps_elem_face;
172 0 : _an = _cs / fi.neighborPtr()->n_sides() * rho_face / eps_elem_face;
173 : }
174 : }
175 : else
176 : {
177 : const auto [interp_coeffs, advected] =
178 187127729 : interpCoeffsAndAdvected(*_rho_u, advected_face_arg, state);
179 :
180 187127729 : const auto elem_arg = elemArg();
181 187127729 : const auto neighbor_arg = neighborArg();
182 :
183 187127729 : const auto rho_elem = _rho(elem_arg, state), rho_neighbor = _rho(neighbor_arg, state);
184 187127729 : const auto eps_elem = epsilon()(elem_arg, state),
185 187127729 : eps_neighbor = epsilon()(neighbor_arg, state);
186 187127729 : const auto var_elem = advected.first / rho_elem * eps_elem,
187 187127729 : var_neighbor = advected.second / rho_neighbor * eps_neighbor;
188 :
189 187127729 : _ae = vdotn * rho_elem / eps_elem * interp_coeffs.first;
190 : // Minus sign because we apply a minus sign to the residual in computeResidual
191 374255458 : _an = -vdotn * rho_neighbor / eps_neighbor * interp_coeffs.second;
192 :
193 187127729 : _elem_residual = _ae * var_elem - _an * var_neighbor;
194 187127729 : _neighbor_residual = -_elem_residual;
195 :
196 187127729 : if (_approximate_as)
197 : {
198 12631840 : _ae = _cs / fi.elem().n_sides() * rho_elem / eps_elem;
199 18947760 : _an = _cs / fi.neighborPtr()->n_sides() * rho_neighbor / eps_neighbor;
200 : }
201 : }
202 : }
203 :
204 190408499 : _ae *= fi.faceArea() * fi.faceCoord();
205 190408499 : _an *= fi.faceArea() * fi.faceCoord();
206 190408499 : _elem_residual *= fi.faceArea() * fi.faceCoord();
207 190408499 : _neighbor_residual *= fi.faceArea() * fi.faceCoord();
208 190408499 : }
209 :
210 : void
211 54298192 : INSFVMomentumAdvection::computeResidual(const FaceInfo & fi)
212 : {
213 54298192 : if (skipForBoundary(fi))
214 : return;
215 :
216 50911248 : computeResidualsAndAData(fi);
217 :
218 50911248 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
219 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
220 : {
221 : // residual contribution of this kernel to the elem element
222 50911098 : prepareVectorTag(_assembly, _var.number());
223 50911098 : _local_re(0) = MetaPhysicL::raw_value(_elem_residual);
224 50911098 : accumulateTaggedLocalResidual();
225 : }
226 :
227 50911248 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
228 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
229 : {
230 : // residual contribution of this kernel to the neighbor element
231 50118768 : prepareVectorTagNeighbor(_assembly, _var.number());
232 50118768 : _local_re(0) = MetaPhysicL::raw_value(_neighbor_residual);
233 50118768 : accumulateTaggedLocalResidual();
234 : }
235 : }
236 :
237 : void
238 54883210 : INSFVMomentumAdvection::computeJacobian(const FaceInfo & fi)
239 : {
240 54883210 : if (skipForBoundary(fi))
241 : return;
242 :
243 49773844 : computeResidualsAndAData(fi);
244 :
245 49773844 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
246 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
247 : {
248 : mooseAssert(_var.dofIndices().size() == 1, "We're currently built to use CONSTANT MONOMIALS");
249 :
250 49773744 : addResidualsAndJacobian(_assembly,
251 99547488 : std::array<ADReal, 1>{{_elem_residual}},
252 : _var.dofIndices(),
253 49773744 : _var.scalingFactor());
254 : }
255 :
256 49773844 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
257 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
258 : {
259 : mooseAssert((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) ==
260 : (_var.dofIndices().size() == 0),
261 : "If the variable is only defined on the neighbor hand side of the face, then that "
262 : "means it should have no dof indices on the elem element. Conversely if "
263 : "the variable is defined on both sides of the face, then it should have a non-zero "
264 : "number of degrees of freedom on the elem element");
265 :
266 : mooseAssert(_var.dofIndicesNeighbor().size() == 1,
267 : "We're currently built to use CONSTANT MONOMIALS");
268 :
269 48692456 : addResidualsAndJacobian(_assembly,
270 97384912 : std::array<ADReal, 1>{{_neighbor_residual}},
271 : _var.dofIndicesNeighbor(),
272 48692456 : _var.scalingFactor());
273 : }
274 : }
275 :
276 : void
277 95771881 : INSFVMomentumAdvection::gatherRCData(const FaceInfo & fi)
278 : {
279 95771881 : if (skipForBoundary(fi))
280 : return;
281 :
282 89865391 : const auto saved_velocity_interp_method = _velocity_interp_method;
283 89865391 : _velocity_interp_method = Moose::FV::InterpMethod::Average;
284 : // Fill-in the coefficients _ae and _an
285 89865391 : computeResidualsAndAData(fi);
286 89865391 : _velocity_interp_method = saved_velocity_interp_method;
287 :
288 89865391 : if (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
289 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
290 89865141 : _rc_uo.addToA(&fi.elem(), _index, _ae);
291 89865391 : if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
292 : _face_type == FaceInfo::VarFaceNeighbors::BOTH)
293 176950954 : _rc_uo.addToA(fi.neighborPtr(), _index, _an);
294 : }
|
