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 40220174 : loopOverElemFaceInfo(const Elem & elem,
42 : const MooseMesh & mesh,
43 : const SubProblem & subproblem,
44 : ActionFunctor & act)
45 : {
46 40220174 : const auto coord_type = subproblem.getCoordSystem(elem.subdomain_id());
47 43084622 : loopOverElemFaceInfo(elem,
48 : mesh,
49 : act,
50 : coord_type,
51 2864448 : coord_type == Moose::COORD_RZ ? subproblem.getAxisymmetricRadialCoord()
52 : : libMesh::invalid_uint);
53 40220174 : }
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 4550854 : 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 13652562 : [this, &fi, elem_to_extrapolate_from]() -> std::pair<const Elem *, bool>
80 : {
81 4550854 : if (elem_to_extrapolate_from)
82 4073378 : 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 4550854 : }();
93 :
94 4550854 : if (two_term_expansion && isFullyDevelopedFlowFace(fi))
95 : {
96 : const Point vector_to_face = elem_to_extrapolate_from_is_fi_elem
97 3595254 : ? (fi.faceCentroid() - fi.elemCentroid())
98 : : (fi.faceCentroid() - fi.neighborCentroid());
99 3595254 : boundary_value = uncorrectedAdGradSln(fi, time, correct_skewness) * vector_to_face +
100 7190508 : getElemValue(elem_to_extrapolate_from, time);
101 : }
102 : else
103 1911200 : boundary_value = INSFVVariable::getExtrapolatedBoundaryFaceValue(
104 955600 : fi, two_term_expansion, correct_skewness, elem_to_extrapolate_from, time);
105 :
106 4550854 : return boundary_value;
107 : }
108 :
109 : VectorValue<ADReal>
110 105185424 : INSFVVelocityVariable::uncorrectedAdGradSln(const FaceInfo & fi,
111 : const StateArg & time,
112 : const bool correct_skewness) const
113 : {
114 : VectorValue<ADReal> unc_grad;
115 105185424 : if (_two_term_boundary_expansion && isFullyDevelopedFlowFace(fi))
116 : {
117 3639426 : const auto & cell_gradient = adGradSln(&fi.elem(), time, correct_skewness);
118 : const auto & normal = fi.normal();
119 10918278 : unc_grad = cell_gradient - (cell_gradient * normal) * normal;
120 : }
121 : else
122 203091996 : unc_grad = INSFVVariable::uncorrectedAdGradSln(fi, time, correct_skewness);
123 :
124 105185424 : return unc_grad;
125 : }
126 :
127 : const VectorValue<ADReal> &
128 212658440 : INSFVVelocityVariable::adGradSln(const Elem * const elem,
129 : const StateArg & time,
130 : bool correct_skewness) const
131 : {
132 212658440 : VectorValue<ADReal> * value_pointer = &_temp_cell_gradient;
133 :
134 : // We ensure that no caching takes place when we compute skewness-corrected
135 : // quantities.
136 212658440 : if (_cache_cell_gradients && !correct_skewness && time.state == 0)
137 : {
138 : auto it = _elem_to_grad.find(elem);
139 :
140 212534420 : if (it != _elem_to_grad.end())
141 172438266 : return it->second;
142 : }
143 :
144 40220174 : 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 40220174 : bool volume_set = false;
156 40220174 : 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 40220174 : 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 161737844 : 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 161737844 : if (isExtrapolatedBoundaryFace(*fi, &functor_elem, time))
233 : {
234 2610391 : if (_two_term_boundary_expansion)
235 : {
236 1723010 : const bool fdf_face = isFullyDevelopedFlowFace(*fi);
237 1723010 : ebf_faces.push_back(std::make_pair(fi, fdf_face));
238 :
239 : // eqn. 2
240 3446020 : ebf_grad_coeffs.push_back(-1. * (elem_has_info
241 1723010 : ? (fi->faceCentroid() - fi->elemCentroid())
242 : : (fi->faceCentroid() - fi->neighborCentroid())));
243 1723010 : ebf_b.push_back(&elem_value);
244 :
245 : // eqn. 1
246 1723010 : grad_ebf_coeffs.push_back(-surface_vector);
247 :
248 : // eqn. 3
249 1723010 : if (fdf_face)
250 : {
251 : // Will be nice in C++17 we'll get a returned reference from this method
252 487806 : fdf_grad_centroid_coeffs.emplace_back();
253 487806 : auto & current_coeffs = fdf_grad_centroid_coeffs.back();
254 487806 : const auto normal = fi->normal();
255 1951224 : for (const auto i : make_range(lm_dim))
256 5853672 : for (const auto j : make_range(lm_dim))
257 : {
258 : auto & current_coeff = current_coeffs(i, j);
259 4390254 : current_coeff = normal(i) * normal(j);
260 4390254 : if (i == j)
261 1463418 : 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 887381 : grad_b += surface_vector * elem_value;
271 : }
272 159127453 : else if (isInternalFace(*fi))
273 : grad_b +=
274 307138882 : surface_vector * (*this)(Moose::FaceArg({fi,
275 : Moose::FV::LimiterType::CentralDifference,
276 : true,
277 : correct_skewness,
278 : nullptr,
279 : nullptr}),
280 153569441 : time);
281 : else
282 : {
283 : mooseAssert(isDirichletBoundaryFace(*fi, &functor_elem, time),
284 : "We've run out of face types");
285 5558012 : grad_b += surface_vector * getDirichletBoundaryFaceValue(*fi, &functor_elem, time);
286 : }
287 :
288 161737844 : 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 40220174 : if (elem_has_info)
293 : {
294 2794841 : coordTransformFactor(
295 2794841 : this->_subproblem, functor_elem.subdomain_id(), fi->elemCentroid(), coord);
296 2794841 : volume = fi->elemVolume() * coord;
297 : }
298 : else
299 : {
300 37425333 : coordTransformFactor(
301 37425333 : this->_subproblem, neighbor->subdomain_id(), fi->neighborCentroid(), coord);
302 37425333 : volume = fi->neighborVolume() * coord;
303 : }
304 :
305 40220174 : volume_set = true;
306 : }
307 201958018 : };
308 :
309 40220174 : Moose::FV::loopOverElemFaceInfo(*elem, this->_mesh, this->_subproblem, action_functor);
310 :
311 : mooseAssert(volume_set && volume > 0, "We should have set the volume");
312 40220174 : grad_b /= volume;
313 :
314 40220174 : const auto coord_system = this->_subproblem.getCoordSystem(elem->subdomain_id());
315 40220174 : if (coord_system == Moose::CoordinateSystemType::COORD_RZ)
316 : {
317 2864448 : const auto r_coord = this->_subproblem.getAxisymmetricRadialCoord();
318 5728896 : 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 40220174 : const auto num_ebfs = static_cast<unsigned int>(ebf_faces.size());
331 :
332 : // test for simple case
333 40220174 : 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 1648239 : lm_dim + num_ebfs + lm_dim * static_cast<unsigned int>(fdf_grad_centroid_coeffs.size());
340 1648239 : DenseVector<ADReal> x(sys_dim), b(sys_dim);
341 1648239 : DenseMatrix<ADReal> A(sys_dim, sys_dim);
342 :
343 : // eqn. 1
344 6592956 : for (const auto lm_dim_index : make_range(lm_dim))
345 : {
346 : // LHS term 1 coeffs
347 4944717 : A(lm_dim_index, lm_dim_index) = 1;
348 :
349 : // LHS term 2 coeffs
350 10113747 : for (const auto ebf_index : make_range(num_ebfs))
351 5169030 : A(lm_dim_index, lm_dim + ebf_index) = grad_ebf_coeffs[ebf_index](lm_dim_index) / volume;
352 :
353 : // RHS
354 4944717 : b(lm_dim_index) = grad_b(lm_dim_index);
355 : }
356 :
357 : unsigned int num_fdf_faces = 0;
358 :
359 : // eqn. 2
360 3371249 : for (const auto ebf_index : make_range(num_ebfs))
361 : {
362 : // LHS term 1 coeffs
363 1723010 : A(lm_dim + ebf_index, lm_dim + ebf_index) = 1;
364 :
365 1723010 : const bool fdf_face = ebf_faces[ebf_index].second;
366 : const unsigned int starting_j_index =
367 1723010 : fdf_face ? lm_dim + num_ebfs + num_fdf_faces * lm_dim : 0;
368 :
369 1723010 : num_fdf_faces += fdf_face;
370 :
371 : // LHS term 2 coeffs
372 6892040 : for (const auto lm_dim_index : make_range(lm_dim))
373 5169030 : A(lm_dim + ebf_index, starting_j_index + lm_dim_index) =
374 5169030 : ebf_grad_coeffs[ebf_index](lm_dim_index);
375 :
376 : // RHS
377 1723010 : 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 2136045 : for (const auto fdf_face_index : make_range(num_fdf_faces))
386 : {
387 487806 : const auto starting_i_index = lm_dim + num_ebfs + fdf_face_index * lm_dim;
388 :
389 1951224 : for (const auto lm_dim_i_index : make_range(lm_dim))
390 : {
391 1463418 : auto i_index = starting_i_index + lm_dim_i_index;
392 1463418 : A(i_index, i_index) = 1;
393 :
394 5853672 : 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 4390254 : fdf_grad_centroid_coeffs[fdf_face_index](lm_dim_i_index, lm_dim_j_index);
398 : }
399 : }
400 :
401 1648239 : A.lu_solve(b, x);
402 6586212 : for (const auto lm_dim_index : make_range(lm_dim))
403 4939659 : grad(lm_dim_index) = x(lm_dim_index);
404 1648239 : }
405 :
406 40218488 : 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 40094468 : return pr.first->second;
411 : }
412 : else
413 : return grad;
414 40225232 : }
415 1686 : catch (std::exception & e)
416 : {
417 : // Don't ignore any *real* errors; we only handle matrix
418 : // singularity here
419 1686 : 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 1686 : const_cast<INSFVVelocityVariable *>(this)->_two_term_boundary_expansion = false;
427 1686 : 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 1686 : const_cast<INSFVVelocityVariable *>(this)->_two_term_boundary_expansion = true;
432 :
433 : return grad;
434 1686 : }
435 40220174 : }
|
