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 "INSFVVelocityVariable.h"
11 : #include "FaceInfo.h"
12 : #include "ADReal.h"
13 : #include "MathFVUtils.h"
14 : #include "FVUtils.h"
15 : #include "MooseError.h"
16 : #include "SubProblem.h"
17 : #include "INSFVNoSlipWallBC.h"
18 : #include "INSFVAttributes.h"
19 : #include "SystemBase.h"
20 : #include "FVDirichletBCBase.h"
21 : #include "Assembly.h"
22 :
23 : #include "libmesh/elem.h"
24 : #include "libmesh/vector_value.h"
25 : #include "libmesh/tensor_value.h"
26 : #include "libmesh/dense_vector.h"
27 : #include "libmesh/dense_matrix.h"
28 : #include "libmesh/libmesh_common.h"
29 :
30 : #include <vector>
31 : #include <utility>
32 :
33 : namespace Moose
34 : {
35 : namespace FV
36 : {
37 : template <typename ActionFunctor>
38 : void
39 55873265 : loopOverElemFaceInfo(const Elem & elem,
40 : const MooseMesh & mesh,
41 : const SubProblem & subproblem,
42 : ActionFunctor & act)
43 : {
44 55873265 : const auto coord_type = subproblem.getCoordSystem(elem.subdomain_id());
45 60199371 : loopOverElemFaceInfo(elem,
46 : mesh,
47 : act,
48 : coord_type,
49 4326106 : coord_type == Moose::COORD_RZ ? subproblem.getAxisymmetricRadialCoord()
50 : : libMesh::invalid_uint);
51 55873265 : }
52 : }
53 : }
54 :
55 : registerMooseObject("NavierStokesApp", INSFVVelocityVariable);
56 :
57 : InputParameters
58 32780 : INSFVVelocityVariable::validParams()
59 : {
60 32780 : return INSFVVariable::validParams();
61 : }
62 :
63 17393 : INSFVVelocityVariable::INSFVVelocityVariable(const InputParameters & params) : INSFVVariable(params)
64 : {
65 17393 : }
66 :
67 : ADReal
68 6183569 : INSFVVelocityVariable::getExtrapolatedBoundaryFaceValue(const FaceInfo & fi,
69 : const bool two_term_expansion,
70 : const bool correct_skewness,
71 : const Elem * elem_to_extrapolate_from,
72 : const StateArg & time) const
73 : {
74 : ADReal boundary_value;
75 : bool elem_to_extrapolate_from_is_fi_elem;
76 : std::tie(elem_to_extrapolate_from, elem_to_extrapolate_from_is_fi_elem) =
77 18550707 : [this, &fi, elem_to_extrapolate_from]() -> std::pair<const Elem *, bool>
78 : {
79 6183569 : if (elem_to_extrapolate_from)
80 5637471 : return {elem_to_extrapolate_from, elem_to_extrapolate_from == &fi.elem()};
81 : else
82 : {
83 : const auto [elem_guaranteed_to_have_dofs,
84 : other_elem,
85 : elem_guaranteed_to_have_dofs_is_fi_elem] =
86 546098 : Moose::FV::determineElemOneAndTwo(fi, *this);
87 : libmesh_ignore(other_elem);
88 546098 : return {elem_guaranteed_to_have_dofs, elem_guaranteed_to_have_dofs_is_fi_elem};
89 : }
90 6183569 : }();
91 :
92 6183569 : if (two_term_expansion && isFullyDevelopedFlowFace(fi))
93 : {
94 : const Point vector_to_face = elem_to_extrapolate_from_is_fi_elem
95 4892238 : ? (fi.faceCentroid() - fi.elemCentroid())
96 : : (fi.faceCentroid() - fi.neighborCentroid());
97 4892238 : boundary_value = uncorrectedAdGradSln(fi, time, correct_skewness) * vector_to_face +
98 9784476 : getElemValue(elem_to_extrapolate_from, time);
99 : }
100 : else
101 2582662 : boundary_value = INSFVVariable::getExtrapolatedBoundaryFaceValue(
102 1291331 : fi, two_term_expansion, correct_skewness, elem_to_extrapolate_from, time);
103 :
104 6183569 : return boundary_value;
105 : }
106 :
107 : VectorValue<ADReal>
108 141941522 : INSFVVelocityVariable::uncorrectedAdGradSln(const FaceInfo & fi,
109 : const StateArg & time,
110 : const bool correct_skewness) const
111 : {
112 : VectorValue<ADReal> unc_grad;
113 141941522 : if (_two_term_boundary_expansion && isFullyDevelopedFlowFace(fi))
114 : {
115 4958664 : const auto & cell_gradient = adGradSln(&fi.elem(), time, correct_skewness);
116 : const auto & normal = fi.normal();
117 14875992 : unc_grad = cell_gradient - (cell_gradient * normal) * normal;
118 : }
119 : else
120 273965716 : unc_grad = INSFVVariable::uncorrectedAdGradSln(fi, time, correct_skewness);
121 :
122 141941522 : return unc_grad;
123 : }
124 :
125 : const VectorValue<ADReal> &
126 285783177 : INSFVVelocityVariable::adGradSln(const Elem * const elem,
127 : const StateArg & time,
128 : bool correct_skewness) const
129 : {
130 285783177 : VectorValue<ADReal> * value_pointer = &_temp_cell_gradient;
131 :
132 : // We ensure that no caching takes place when we compute skewness-corrected
133 : // quantities.
134 285783177 : if (_cache_cell_gradients && !correct_skewness && time.state == 0)
135 : {
136 : auto it = _elem_to_grad.find(elem);
137 :
138 285597147 : if (it != _elem_to_grad.end())
139 229909912 : return it->second;
140 : }
141 :
142 55873265 : ADReal elem_value = getElemValue(elem, time);
143 :
144 : // We'll save off the extrapolated boundary faces (ebf) for later assignment to the cache (these
145 : // are the keys). The boolean in the pair will denote whether the ebf face is a fully developed
146 : // flow (e.g. fdf) face
147 : std::vector<std::pair<const FaceInfo *, bool>> ebf_faces;
148 :
149 : try
150 : {
151 : VectorValue<ADReal> & grad = *value_pointer;
152 :
153 55873265 : bool volume_set = false;
154 55873265 : Real volume = 0;
155 :
156 : // If we are performing a two term Taylor expansion for extrapolated boundary faces (faces on
157 : // boundaries that do not have associated Dirichlet conditions), then the element gradient
158 : // depends on the boundary face value and the boundary face value depends on the element
159 : // gradient, so we have a system of equations to solve. Here is the system:
160 : //
161 : // \nabla \phi_C - \frac{1}{V} \sum_{ebf} \phi_{ebf} \vec{S_f} =
162 : // \frac{1}{V} \sum_{of} \phi_{of} \vec{S_f} eqn.
163 : // 1
164 : //
165 : // \phi_{ebf} - \vec{d_{Cf}} \cdot \nabla \phi_C = \phi_C eqn.
166 : // 2
167 : //
168 : // where $C$ refers to the cell centroid, $ebf$ refers to an extrapolated boundary face, $of$
169 : // refers to "other faces", e.g. non-ebf faces, and $f$ is a general face. $d_{Cf}$ is the
170 : // vector drawn from the element centroid to the face centroid, and $\vec{S_f}$ is the surface
171 : // vector, e.g. the face area times the outward facing normal
172 : //
173 : // NOTE: On fully developed flow boundaries, we modify our equation set slightly. In equation
174 : // 2,
175 : // $\nabla \phi_C$ is replaced with $\nabla \phi_{ebf,fdf}$ where $fdf$ denotes fully
176 : // developed flow. Moreover, we introduce a third equation:
177 : //
178 : // \nabla \phi_{ebf,fdf} - \nabla \phi_C + (\nabla \phi_C \cdot \hat{n}) \hat{n} = 0 eqn.
179 : // 3
180 : //
181 : // These modifications correspond to Moukalled's equations 15.140 and 15.141, but with
182 : // $\hat{e_b}$ replaced with $\hat{n}$ because we believe the equation as written doesn't
183 : // reflect the intent of the text, which is to guarantee a zero normal gradient in the
184 : // direction of the surface normal
185 :
186 : // ebf eqns: element gradient coefficients, e.g. eqn. 2, LHS term 2 coefficient. *Note* that
187 : // each element of the std::vector could correspond to a cell centroid gradient or to a face
188 : // gradient computed on a fully developed flow face
189 : std::vector<VectorValue<Real>> ebf_grad_coeffs;
190 : // ebf eqns: rhs b values. These will actually correspond to the elem_value so we can use a
191 : // pointer and avoid copying. This is the RHS of eqn. 2
192 : std::vector<const ADReal *> ebf_b;
193 :
194 : // elem grad eqns: ebf coefficients, e.g. eqn. 1, LHS term 2 coefficients
195 : std::vector<VectorValue<Real>> grad_ebf_coeffs;
196 : // elem grad eqns: rhs b value, e.g. eqn. 1 RHS
197 55873265 : VectorValue<ADReal> grad_b = 0;
198 :
199 : // eqn. 3 coefficients for cell centroid gradient, e.g. the coefficients that fall out of term
200 : // 2 on the LHS of eqn. 3
201 : std::vector<TensorValue<Real>> fdf_grad_centroid_coeffs;
202 :
203 : const unsigned int lm_dim = LIBMESH_DIM;
204 :
205 224940674 : auto action_functor = [&volume_set,
206 : &volume,
207 : &elem_value,
208 : #ifndef NDEBUG
209 : &elem,
210 : #endif
211 : &ebf_faces,
212 : &ebf_grad_coeffs,
213 : &ebf_b,
214 : &grad_ebf_coeffs,
215 : &grad_b,
216 : &fdf_grad_centroid_coeffs,
217 : correct_skewness,
218 : &time,
219 : this](const Elem & functor_elem,
220 : const Elem * const neighbor,
221 : const FaceInfo * const fi,
222 : const Point & surface_vector,
223 : Real coord,
224 : const bool elem_has_info)
225 : {
226 : mooseAssert(fi, "We need a FaceInfo for this action_functor");
227 : mooseAssert(elem == &functor_elem,
228 : "Just a sanity check that the element being passed in is the one we passed out.");
229 :
230 224940674 : if (isExtrapolatedBoundaryFace(*fi, &functor_elem, time))
231 : {
232 3504906 : if (_two_term_boundary_expansion)
233 : {
234 2319249 : const bool fdf_face = isFullyDevelopedFlowFace(*fi);
235 2319249 : ebf_faces.push_back(std::make_pair(fi, fdf_face));
236 :
237 : // eqn. 2
238 4638498 : ebf_grad_coeffs.push_back(-1. * (elem_has_info
239 2319249 : ? (fi->faceCentroid() - fi->elemCentroid())
240 : : (fi->faceCentroid() - fi->neighborCentroid())));
241 2319249 : ebf_b.push_back(&elem_value);
242 :
243 : // eqn. 1
244 2319249 : grad_ebf_coeffs.push_back(-surface_vector);
245 :
246 : // eqn. 3
247 2319249 : if (fdf_face)
248 : {
249 : // Will be nice in C++17 we'll get a returned reference from this method
250 675001 : fdf_grad_centroid_coeffs.emplace_back();
251 675001 : auto & current_coeffs = fdf_grad_centroid_coeffs.back();
252 675001 : const auto normal = fi->normal();
253 2700004 : for (const auto i : make_range(lm_dim))
254 8100012 : for (const auto j : make_range(lm_dim))
255 : {
256 : auto & current_coeff = current_coeffs(i, j);
257 6075009 : current_coeff = normal(i) * normal(j);
258 6075009 : if (i == j)
259 2025003 : current_coeff -= 1.;
260 : }
261 : }
262 : }
263 : else
264 : // We are doing a one-term expansion for the extrapolated boundary faces, in which case
265 : // we have no eqn. 2 and we have no second term in the LHS of eqn. 1. Instead we apply
266 : // the element centroid value as the face value (one-term expansion) in the RHS of eqn.
267 : // 1
268 1185657 : grad_b += surface_vector * elem_value;
269 : }
270 221435768 : else if (isInternalFace(*fi))
271 : grad_b +=
272 426914834 : surface_vector * (*this)(Moose::FaceArg({fi,
273 : Moose::FV::LimiterType::CentralDifference,
274 : true,
275 : correct_skewness,
276 : nullptr,
277 : nullptr}),
278 213457417 : time);
279 : else
280 : {
281 : mooseAssert(isDirichletBoundaryFace(*fi, &functor_elem, time),
282 : "We've run out of face types");
283 7978351 : grad_b += surface_vector * getDirichletBoundaryFaceValue(*fi, &functor_elem, time);
284 : }
285 :
286 224940674 : if (!volume_set)
287 : {
288 : // We use the FaceInfo volumes because those values have been pre-computed and cached.
289 : // An explicit call to elem->volume() here would incur unnecessary expense
290 55873265 : if (elem_has_info)
291 : {
292 3882673 : coordTransformFactor(
293 3882673 : this->_subproblem, functor_elem.subdomain_id(), fi->elemCentroid(), coord);
294 3882673 : volume = fi->elemVolume() * coord;
295 : }
296 : else
297 : {
298 51990592 : coordTransformFactor(
299 51990592 : this->_subproblem, neighbor->subdomain_id(), fi->neighborCentroid(), coord);
300 51990592 : volume = fi->neighborVolume() * coord;
301 : }
302 :
303 55873265 : volume_set = true;
304 : }
305 280813939 : };
306 :
307 55873265 : Moose::FV::loopOverElemFaceInfo(*elem, this->_mesh, this->_subproblem, action_functor);
308 :
309 : mooseAssert(volume_set && volume > 0, "We should have set the volume");
310 55873265 : grad_b /= volume;
311 :
312 55873265 : const auto coord_system = this->_subproblem.getCoordSystem(elem->subdomain_id());
313 55873265 : if (coord_system == Moose::CoordinateSystemType::COORD_RZ)
314 : {
315 4326106 : const auto r_coord = this->_subproblem.getAxisymmetricRadialCoord();
316 8652212 : grad_b(r_coord) -= elem_value / elem->vertex_average()(r_coord);
317 : }
318 :
319 : mooseAssert(
320 : coord_system != Moose::CoordinateSystemType::COORD_RSPHERICAL,
321 : "We have not yet implemented the correct translation from gradient to divergence for "
322 : "spherical coordinates yet.");
323 :
324 : mooseAssert(
325 : ebf_faces.size() < UINT_MAX,
326 : "You've created a mystical element that has more faces than can be held by unsigned "
327 : "int. I applaud you.");
328 55873265 : const auto num_ebfs = static_cast<unsigned int>(ebf_faces.size());
329 :
330 : // test for simple case
331 55873265 : if (num_ebfs == 0)
332 : grad = grad_b;
333 : else
334 : {
335 : // We have to solve a system
336 : const unsigned int sys_dim =
337 2209910 : lm_dim + num_ebfs + lm_dim * static_cast<unsigned int>(fdf_grad_centroid_coeffs.size());
338 2209910 : DenseVector<ADReal> x(sys_dim), b(sys_dim);
339 2209910 : DenseMatrix<ADReal> A(sys_dim, sys_dim);
340 :
341 : // eqn. 1
342 8839640 : for (const auto lm_dim_index : make_range(lm_dim))
343 : {
344 : // LHS term 1 coeffs
345 6629730 : A(lm_dim_index, lm_dim_index) = 1;
346 :
347 : // LHS term 2 coeffs
348 13587477 : for (const auto ebf_index : make_range(num_ebfs))
349 6957747 : A(lm_dim_index, lm_dim + ebf_index) = grad_ebf_coeffs[ebf_index](lm_dim_index) / volume;
350 :
351 : // RHS
352 6629730 : b(lm_dim_index) = grad_b(lm_dim_index);
353 : }
354 :
355 : unsigned int num_fdf_faces = 0;
356 :
357 : // eqn. 2
358 4529159 : for (const auto ebf_index : make_range(num_ebfs))
359 : {
360 : // LHS term 1 coeffs
361 2319249 : A(lm_dim + ebf_index, lm_dim + ebf_index) = 1;
362 :
363 2319249 : const bool fdf_face = ebf_faces[ebf_index].second;
364 : const unsigned int starting_j_index =
365 2319249 : fdf_face ? lm_dim + num_ebfs + num_fdf_faces * lm_dim : 0;
366 :
367 2319249 : num_fdf_faces += fdf_face;
368 :
369 : // LHS term 2 coeffs
370 9276996 : for (const auto lm_dim_index : make_range(lm_dim))
371 6957747 : A(lm_dim + ebf_index, starting_j_index + lm_dim_index) =
372 6957747 : ebf_grad_coeffs[ebf_index](lm_dim_index);
373 :
374 : // RHS
375 2319249 : b(lm_dim + ebf_index) = *ebf_b[ebf_index];
376 : }
377 :
378 : mooseAssert(num_fdf_faces == fdf_grad_centroid_coeffs.size(),
379 : "Bad math in INSFVVelocityVariable::adGradlnSln(const Elem *). Please contact a "
380 : "MOOSE developer");
381 :
382 : // eqn. 3
383 2884911 : for (const auto fdf_face_index : make_range(num_fdf_faces))
384 : {
385 675001 : const auto starting_i_index = lm_dim + num_ebfs + fdf_face_index * lm_dim;
386 :
387 2700004 : for (const auto lm_dim_i_index : make_range(lm_dim))
388 : {
389 2025003 : auto i_index = starting_i_index + lm_dim_i_index;
390 2025003 : A(i_index, i_index) = 1;
391 :
392 8100012 : for (const auto lm_dim_j_index : make_range(lm_dim))
393 : // j_index = lm_dim_j_index
394 : A(i_index, lm_dim_j_index) =
395 6075009 : fdf_grad_centroid_coeffs[fdf_face_index](lm_dim_i_index, lm_dim_j_index);
396 : }
397 : }
398 :
399 2209910 : A.lu_solve(b, x);
400 8827960 : for (const auto lm_dim_index : make_range(lm_dim))
401 6620970 : grad(lm_dim_index) = x(lm_dim_index);
402 2209910 : }
403 :
404 55870345 : if (_cache_cell_gradients && !correct_skewness)
405 : {
406 : auto pr = _elem_to_grad.emplace(elem, std::move(grad));
407 : mooseAssert(pr.second, "Insertion should have just happened.");
408 55684315 : return pr.first->second;
409 : }
410 : else
411 : return grad;
412 55882025 : }
413 2920 : catch (libMesh::LogicError &)
414 : {
415 : // Retry without two-term
416 : mooseAssert(_two_term_boundary_expansion,
417 : "I believe we should only get singular systems when two-term boundary expansion is "
418 : "being used");
419 2920 : const_cast<INSFVVelocityVariable *>(this)->_two_term_boundary_expansion = false;
420 2920 : const auto & grad = adGradSln(elem, time, correct_skewness);
421 :
422 : // Two term boundary expansion should only fail at domain corners. We want to keep trying it
423 : // at other boundary locations
424 2920 : const_cast<INSFVVelocityVariable *>(this)->_two_term_boundary_expansion = true;
425 :
426 : return grad;
427 2920 : }
428 55873265 : }
|
