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