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 : #pragma once
11 :
12 : #include "MooseFunctor.h"
13 : #include "MathFVUtils.h"
14 : #include "FVUtils.h"
15 : #include "MooseMeshUtils.h"
16 : #include "VectorComponentFunctor.h"
17 : #include "ArrayComponentFunctor.h"
18 : #include "libmesh/elem.h"
19 :
20 : namespace Moose
21 : {
22 : namespace FV
23 : {
24 : /**
25 : * Compute a cell gradient using the method of Green-Gauss
26 : * @param elem_arg An element argument specifying the current element and whether to perform skew
27 : * corrections
28 : * @param state_arg A state argument that indicates what temporal / solution iteration data to use
29 : * when evaluating the provided functor
30 : * @param functor The functor that will provide information such as cell and face value evaluations
31 : * necessary to construct the cell gradient
32 : * @param two_term_boundary_expansion Whether to perform a two-term expansion to compute
33 : * extrapolated boundary face values. If this is true, then an implicit system has to be solved. If
34 : * false, then the cell center value will be used as the extrapolated boundary face value
35 : * @param mesh The mesh on which we are computing the gradient
36 : * @return The computed cell gradient
37 : */
38 : template <typename T,
39 : typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
40 : libMesh::VectorValue<T>
41 21011348 : greenGaussGradient(const ElemArg & elem_arg,
42 : const StateArg & state_arg,
43 : const FunctorBase<T> & functor,
44 : const bool two_term_boundary_expansion,
45 : const MooseMesh & mesh,
46 : const bool force_green_gauss = false)
47 : {
48 : mooseAssert(elem_arg.elem, "This should be non-null");
49 21011348 : const auto coord_type = mesh.getCoordSystem(elem_arg.elem->subdomain_id());
50 21011348 : const auto rz_radial_coord = mesh.getAxisymmetricRadialCoord();
51 :
52 21011348 : const T elem_value = functor(elem_arg, state_arg);
53 :
54 : // We'll count the number of extrapolated boundary faces (ebfs)
55 21011348 : unsigned int num_ebfs = 0;
56 :
57 : // Gradient to be returned
58 21011348 : VectorValue<T> grad;
59 :
60 21011348 : if (!elem_arg.correct_skewness || force_green_gauss) // Do Green-Gauss
61 : {
62 : try
63 : {
64 19480413 : bool volume_set = false;
65 19480413 : Real volume = 0;
66 :
67 : // If we are performing a two term Taylor expansion for extrapolated boundary faces (faces on
68 : // boundaries that do not have associated Dirichlet conditions), then the element gradient
69 : // depends on the boundary face value and the boundary face value depends on the element
70 : // gradient, so we have a system of equations to solve. Here is the system:
71 : //
72 : // \nabla \phi_C - \frac{1}{V} \sum_{ebf} \phi_{ebf} \vec{S_f} =
73 : // \frac{1}{V} \sum_{of} \phi_{of} \vec{S_f} eqn. 1
74 : //
75 : // \phi_{ebf} - \vec{d_{Cf}} \cdot \nabla \phi_C = \phi_C eqn. 2
76 : //
77 : // where $C$ refers to the cell centroid, $ebf$ refers to an extrapolated boundary face, $of$
78 : // refers to "other faces", e.g. non-ebf faces, and $f$ is a general face. $d_{Cf}$ is the
79 : // vector drawn from the element centroid to the face centroid, and $\vec{S_f}$ is the surface
80 : // vector, e.g. the face area times the outward facing normal
81 :
82 : // ebf eqns: element gradient coefficients, e.g. eqn. 2, LHS term 2 coefficient
83 19480413 : std::vector<libMesh::VectorValue<Real>> ebf_grad_coeffs;
84 : // ebf eqns: rhs b values. These will actually correspond to the elem_value so we can use a
85 : // pointer and avoid copying. This is the RHS of eqn. 2
86 19480413 : std::vector<const T *> ebf_b;
87 :
88 : // elem grad eqns: ebf coefficients, e.g. eqn. 1, LHS term 2 coefficients
89 19480413 : std::vector<libMesh::VectorValue<Real>> grad_ebf_coeffs;
90 : // elem grad eqns: rhs b value, e.g. eqn. 1 RHS
91 19480413 : libMesh::VectorValue<T> grad_b = 0;
92 :
93 78996187 : auto action_functor = [&volume_set,
94 : &volume,
95 : &elem_value,
96 : &elem_arg,
97 : &num_ebfs,
98 : &ebf_grad_coeffs,
99 : &ebf_b,
100 : &grad_ebf_coeffs,
101 : &grad_b,
102 : &state_arg,
103 : &functor,
104 : two_term_boundary_expansion,
105 : coord_type,
106 : rz_radial_coord](const Elem & libmesh_dbg_var(functor_elem),
107 : const Elem *,
108 : const FaceInfo * const fi,
109 : const Point & surface_vector,
110 : Real coord,
111 : const bool elem_has_info)
112 : {
113 : mooseAssert(fi, "We need a FaceInfo for this action_functor");
114 : mooseAssert(
115 : elem_arg.elem == &functor_elem,
116 : "Just a sanity check that the element being passed in is the one we passed out.");
117 :
118 73246005 : if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
119 : {
120 2192245 : if (two_term_boundary_expansion)
121 : {
122 474781 : num_ebfs += 1;
123 :
124 : // eqn. 2
125 474781 : ebf_grad_coeffs.push_back(-1. * (elem_has_info
126 474781 : ? (fi->faceCentroid() - fi->elemCentroid())
127 19511 : : (fi->faceCentroid() - fi->neighborCentroid())));
128 474781 : ebf_b.push_back(&elem_value);
129 :
130 : // eqn. 1
131 474781 : grad_ebf_coeffs.push_back(-surface_vector);
132 : }
133 : else
134 : // We are doing a one-term expansion for the extrapolated boundary faces, in which case
135 : // we have no eqn. 2 and we have no second term in the LHS of eqn. 1. Instead we apply
136 : // the element centroid value as the face value (one-term expansion) in the RHS of eqn.
137 : // 1
138 1717464 : grad_b += surface_vector * elem_value;
139 : }
140 : else
141 71053760 : grad_b +=
142 880 : surface_vector * functor(Moose::FaceArg{fi,
143 : Moose::FV::LimiterType::CentralDifference,
144 : true,
145 71053760 : elem_arg.correct_skewness,
146 71053760 : elem_arg.elem,
147 : nullptr},
148 : state_arg);
149 :
150 73246005 : if (!volume_set)
151 : {
152 : // We use the FaceInfo volumes because those values have been pre-computed and cached.
153 : // An explicit call to elem->volume() here would incur unnecessary expense
154 19480413 : if (elem_has_info)
155 : {
156 1380515 : MooseMeshUtils::coordTransformFactor(
157 : fi->elemCentroid(), coord, coord_type, rz_radial_coord);
158 1380515 : volume = fi->elemVolume() * coord;
159 : }
160 : else
161 : {
162 18099898 : MooseMeshUtils::coordTransformFactor(
163 : fi->neighborCentroid(), coord, coord_type, rz_radial_coord);
164 18099898 : volume = fi->neighborVolume() * coord;
165 : }
166 :
167 19480413 : volume_set = true;
168 : }
169 : };
170 :
171 19480413 : Moose::FV::loopOverElemFaceInfo(
172 19480413 : *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
173 :
174 : mooseAssert(volume_set && volume > 0, "We should have set the volume");
175 19480413 : grad_b /= volume;
176 :
177 19480413 : if (coord_type == Moose::CoordinateSystemType::COORD_RZ)
178 : {
179 : mooseAssert(rz_radial_coord != libMesh::invalid_uint, "rz_radial_coord must be set");
180 32365 : grad_b(rz_radial_coord) -= elem_value / elem_arg.elem->vertex_average()(rz_radial_coord);
181 : }
182 :
183 : mooseAssert(
184 : coord_type != Moose::CoordinateSystemType::COORD_RSPHERICAL,
185 : "We have not yet implemented the correct translation from gradient to divergence for "
186 : "spherical coordinates yet.");
187 :
188 : // test for simple case
189 19480413 : if (num_ebfs == 0)
190 19021454 : grad = grad_b;
191 : else
192 : {
193 : // We have to solve a system
194 458959 : const unsigned int sys_dim = Moose::dim + num_ebfs;
195 458959 : DenseVector<T> x(sys_dim), b(sys_dim);
196 458959 : DenseMatrix<T> A(sys_dim, sys_dim);
197 :
198 : // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
199 :
200 : // eqn. 1
201 1835836 : for (const auto i : make_range(Moose::dim))
202 : {
203 : // LHS term 1 coeffs
204 1376877 : A(i, i) = 1;
205 :
206 : // LHS term 2 coeffs
207 2801220 : for (const auto j : make_range(num_ebfs))
208 1424343 : A(i, Moose::dim + j) = grad_ebf_coeffs[j](i) / volume;
209 :
210 : // RHS
211 1376877 : b(i) = grad_b(i);
212 : }
213 :
214 : // eqn. 2
215 933740 : for (const auto j : make_range(num_ebfs))
216 : {
217 : // LHS term 1 coeffs
218 474781 : A(Moose::dim + j, Moose::dim + j) = 1;
219 :
220 : // LHS term 2 coeffs
221 1899124 : for (const auto i : make_range(unsigned(Moose::dim)))
222 1424343 : A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
223 :
224 : // RHS
225 474781 : b(Moose::dim + j) = *ebf_b[j];
226 : }
227 :
228 458959 : A.lu_solve(b, x);
229 : // libMesh is generous about what it considers nonsingular. Let's check a little more
230 : // strictly
231 456805 : if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
232 364 : throw libMesh::LogicError("Matrix A is singular!");
233 1825764 : for (const auto i : make_range(Moose::dim))
234 1369323 : grad(i) = x(i);
235 463995 : }
236 19487967 : }
237 2518 : catch (libMesh::LogicError & e)
238 : {
239 : // Retry without two-term
240 2518 : if (!two_term_boundary_expansion)
241 0 : mooseError(
242 : "I believe we should only get singular systems when two-term boundary expansion is "
243 : "being used. The error thrown during the computation of the gradient: ",
244 0 : e.what());
245 :
246 2518 : grad = greenGaussGradient(elem_arg, state_arg, functor, false, mesh, true);
247 : }
248 19480413 : }
249 : else // Do Least-Squares
250 : {
251 : try
252 : {
253 :
254 : // The least squares method aims to find the gradient at the element centroid by minimizing
255 : // the difference between the estimated and actual value differences across neighboring cells.
256 : // The least squares formulation is:
257 : //
258 : // Minimize J = \sum_{n} [ w_n ((\nabla \phi_C \cdot \delta \vec{x}_n) - \delta \phi_n) ]^2
259 : //
260 : // where:
261 : // - \(\nabla \phi_C\) is the gradient at the element centroid C (unknown).
262 : // - \(\delta \vec{x}_n = \vec{x}_n - \vec{x}_C\) is the vector from the element centroid to
263 : // neighbor \(n\).
264 : // - \(\delta \phi_n = \phi_n - \phi_C\) is the difference in the scalar value between
265 : // neighbor \(n\) and element \(C\).
266 : // - w_n = 1/||\vec{x}_n|| is a vector of weights that is used to handle large aspect ratio
267 : // cells
268 : //
269 : // To handle extrapolated boundary faces (faces on boundaries without Dirichlet conditions),
270 : // we include additional unknowns and equations in the least squares system.
271 : // For ebfs, we may not know the value \(\phi_n\), so we include it as an unknown.
272 : // This results in an augmented system that simultaneously solves for the gradient and the ebf
273 : // values.
274 :
275 : // Get estimated number of faces in a cell
276 1530935 : const auto ptr_range = elem_arg.elem->neighbor_ptr_range();
277 1530935 : const size_t num_faces = std::distance(ptr_range.begin(), ptr_range.end());
278 :
279 : // Lists to store position differences and value differences for least squares computation
280 1530935 : std::vector<Point> delta_x_list; // List of position differences between neighbor centroids
281 : // and element centroid
282 1530935 : delta_x_list.reserve(num_faces);
283 :
284 : std::vector<T>
285 1530935 : delta_phi_list; // List of value differences between neighbor values and element value
286 1530935 : delta_phi_list.reserve(num_faces);
287 :
288 : // Variables to handle extrapolated boundary faces (ebfs) in the least squares method
289 1530935 : std::vector<Point> delta_x_ebf_list; // Position differences for ebfs
290 1530935 : delta_phi_list.reserve(num_faces);
291 :
292 : std::vector<VectorValue<Real>>
293 1530935 : ebf_grad_coeffs; // Coefficients for the gradient in ebf equations
294 1530935 : delta_phi_list.reserve(num_faces);
295 :
296 1530935 : std::vector<const T *> ebf_b; // RHS values for ebf equations (pointers to avoid copying)
297 1530935 : delta_phi_list.reserve(num_faces);
298 :
299 1530935 : unsigned int num_ebfs = 0; // Number of extrapolated boundary faces
300 :
301 : // Action functor to collect data from each face of the element
302 5030215 : auto action_functor = [&elem_value,
303 : &elem_arg,
304 : &num_ebfs,
305 : &ebf_grad_coeffs,
306 : &ebf_b,
307 : &delta_x_list,
308 : &delta_phi_list,
309 : &delta_x_ebf_list,
310 : &state_arg,
311 : &functor,
312 : two_term_boundary_expansion](const Elem & libmesh_dbg_var(loc_elem),
313 : const Elem * loc_neigh,
314 : const FaceInfo * const fi,
315 : const Point & /*surface_vector*/,
316 : Real /*coord*/,
317 : const bool elem_has_info)
318 : {
319 : mooseAssert(fi, "We need a FaceInfo for this action_functor");
320 : mooseAssert(
321 : elem_arg.elem == &loc_elem,
322 : "Just a sanity check that the element being passed in is the one we passed out.");
323 :
324 4592805 : if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
325 : {
326 : // Extrapolated Boundary Face (ebf)
327 17360 : const Point delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
328 0 : : (fi->faceCentroid() - fi->neighborCentroid());
329 17360 : delta_x_list.push_back(delta_x);
330 :
331 17360 : T delta_phi;
332 :
333 17360 : if (two_term_boundary_expansion)
334 : {
335 : // Two-term expansion: include ebf values as unknowns in the augmented system
336 17360 : num_ebfs += 1;
337 :
338 : // Coefficient for the gradient in the ebf equation
339 17360 : ebf_grad_coeffs.push_back(-delta_x);
340 : // RHS value for the ebf equation
341 17360 : ebf_b.push_back(&elem_value);
342 :
343 17360 : delta_phi = -elem_value;
344 17360 : delta_x_ebf_list.push_back(delta_x);
345 : }
346 : else
347 : {
348 : // One-term expansion: assume zero difference across the boundary (\delta \phi = 0)
349 0 : delta_phi = T(0);
350 : }
351 17360 : delta_phi_list.push_back(delta_phi);
352 17360 : }
353 : else
354 : {
355 : // Internal Face or Boundary Face with Dirichlet condition
356 4575445 : Point delta_x;
357 4575445 : T neighbor_value;
358 4575445 : if (functor.isInternalFace(*fi))
359 : {
360 : // Internal face with a neighboring element
361 4494070 : delta_x = loc_neigh->vertex_average() - elem_arg.elem->vertex_average();
362 :
363 4494070 : const ElemArg neighbor_arg = {loc_neigh, elem_arg.correct_skewness};
364 4494070 : neighbor_value = functor(neighbor_arg, state_arg);
365 : }
366 : else
367 : {
368 : // Boundary face with Dirichlet condition
369 81375 : delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
370 0 : : (fi->faceCentroid() - fi->neighborCentroid());
371 81375 : neighbor_value = functor(Moose::FaceArg{fi,
372 : Moose::FV::LimiterType::CentralDifference,
373 : true,
374 81375 : elem_arg.correct_skewness,
375 81375 : elem_arg.elem,
376 : nullptr},
377 : state_arg);
378 : }
379 :
380 4575445 : delta_x_list.push_back(delta_x);
381 4575445 : delta_phi_list.push_back(neighbor_value - elem_value);
382 4575445 : }
383 : };
384 :
385 : // Loop over element faces and apply the action_functor
386 1530935 : Moose::FV::loopOverElemFaceInfo(
387 1530935 : *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
388 :
389 : // Compute Least Squares gradient
390 1530935 : const unsigned int n = delta_x_list.size();
391 : mooseAssert(n >= dim, "Not enough neighbors to perform least squares");
392 :
393 1530935 : DenseMatrix<T> ATA(dim, dim);
394 1530935 : DenseVector<T> ATb(dim);
395 :
396 : // Assemble the normal equations for the least squares method
397 : // ATA = \sum_n (\delta \vec{x}_n^T * \delta \vec{x}_n)
398 : // ATb = \sum_n (\delta \vec{x}_n^T * \delta \phi_n)
399 1530935 : ATA.zero();
400 1530935 : ATb.zero();
401 :
402 6123740 : for (unsigned int i = 0; i < n; ++i)
403 : {
404 4592805 : const Point & delta_x = delta_x_list[i];
405 4592805 : const T & delta_phi = delta_phi_list[i];
406 :
407 18371220 : for (unsigned int d1 = 0; d1 < dim; ++d1)
408 : {
409 13778415 : const Real dx_d1 = delta_x(d1);
410 13778415 : const Real dx_norm = delta_x.norm();
411 13778415 : const Real dx_d1_norm = dx_d1 / dx_norm;
412 :
413 55113660 : for (unsigned int d2 = 0; d2 < dim; ++d2)
414 : {
415 41335245 : const Real dx_d2_norm = delta_x(d2) / dx_norm;
416 41335245 : ATA(d1, d2) += dx_d1_norm * dx_d2_norm;
417 : }
418 :
419 13778415 : ATb(d1) += dx_d1_norm * delta_phi / dx_norm;
420 : }
421 : }
422 :
423 : // Add regularization to ATA to prevent singularity
424 1530935 : T epsilon = 1e-12; // Small regularization parameter
425 6123740 : for (unsigned int d = 0; d < dim; ++d)
426 : {
427 4592805 : ATA(d, d) += epsilon;
428 : }
429 :
430 1530935 : if (num_ebfs == 0)
431 : {
432 : // Solve and store the least squares gradient
433 1513575 : DenseVector<T> grad_ls(dim);
434 1513575 : ATA.lu_solve(ATb, grad_ls);
435 :
436 : // Store the least squares gradient
437 6054300 : for (unsigned int d = 0; d < dim; ++d)
438 : {
439 4540725 : grad(d) = grad_ls(d);
440 : }
441 1513575 : }
442 : else
443 : {
444 : // We have to solve a system of equations
445 17360 : const unsigned int sys_dim = Moose::dim + num_ebfs;
446 17360 : DenseVector<T> x(sys_dim), b(sys_dim);
447 17360 : DenseMatrix<T> A(sys_dim, sys_dim);
448 :
449 : // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
450 :
451 : // eqn. 1: Element gradient equations
452 69440 : for (unsigned int d1 = 0; d1 < dim; ++d1)
453 : {
454 : // LHS term 1 coefficients (gradient components)
455 208320 : for (unsigned int d2 = 0; d2 < dim; ++d2)
456 156240 : A(d1, d2) = ATA(d1, d2);
457 :
458 : // RHS
459 52080 : b(d1) = ATb(d1);
460 : }
461 :
462 : // LHS term 2 coefficients (ebf contributions)
463 34720 : for (const auto i : make_range(num_ebfs))
464 : {
465 17360 : const Point & delta_x = delta_x_ebf_list[i];
466 69440 : for (unsigned int d1 = 0; d1 < dim; ++d1)
467 : {
468 52080 : const Real dx_d1 = delta_x(d1);
469 52080 : A(d1, Moose::dim + i) -= dx_d1 / delta_x.norm() / delta_x.norm();
470 : }
471 : }
472 :
473 : // eqn. 2: Extrapolated boundary face equations
474 34720 : for (const auto j : make_range(num_ebfs))
475 : {
476 : // LHS term 1 coefficients (ebf values)
477 17360 : A(Moose::dim + j, Moose::dim + j) = 1;
478 :
479 : // LHS term 2 coefficients (gradient components)
480 69440 : for (const auto i : make_range(unsigned(Moose::dim)))
481 52080 : A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
482 :
483 : // RHS
484 17360 : b(Moose::dim + j) = *ebf_b[j];
485 : }
486 :
487 : // Solve the system
488 17360 : A.lu_solve(b, x);
489 :
490 : // Check for singularity
491 17360 : if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
492 0 : throw libMesh::LogicError("Matrix A is singular!");
493 :
494 : // Extract the gradient components
495 69440 : for (const auto i : make_range(Moose::dim))
496 52080 : grad(i) = x(i);
497 17360 : }
498 :
499 : mooseAssert(
500 : coord_type != Moose::CoordinateSystemType::COORD_RSPHERICAL,
501 : "We have not yet implemented the correct translation from gradient to divergence for "
502 : "spherical coordinates yet.");
503 1530935 : }
504 0 : catch (libMesh::LogicError & e)
505 : {
506 : // Log warning and default to simple green Gauss
507 0 : mooseWarning(
508 : "Singular matrix encountered in least squares gradient computation. Falling back "
509 : "to Green-Gauss gradient.");
510 :
511 0 : grad = greenGaussGradient(
512 : elem_arg, state_arg, functor, false, mesh, /* force_green_gauss */ true);
513 : }
514 : }
515 :
516 42022696 : return grad;
517 :
518 : // Notes to future developer:
519 : // Note 1:
520 : // For the least squares gradient, the LS matrix could be precomputed and stored for every cell
521 : // I tried doing this on October 2024, but the element lookup for these matrices is too slow
522 : // and seems better to compute weights on the fly.
523 : // Consider building a map from elem_id to these matrices and speed up lookup with octree
524 : // Note 2:
525 : // In the future one would like to have a hybrid gradient scheme, where:
526 : // \nabla \phi = \beta (\nabla \phi)_{LS} + (1 - \beta) (\nabla \phi)_{GG}
527 : // Then optize \beta based on mesh heuristics
528 21011128 : }
529 :
530 : /**
531 : * Compute a face gradient from Green-Gauss cell gradients, with orthogonality correction
532 : * On the boundaries, the boundary element value is used
533 : * @param face_arg A face argument specifying the current faceand whether to perform skew
534 : * corrections
535 : * @param state_arg A state argument that indicates what temporal / solution iteration data to use
536 : * when evaluating the provided functor
537 : * @param functor The functor that will provide information such as cell and face value evaluations
538 : * necessary to construct the cell gradient
539 : * @param two_term_boundary_expansion Whether to perform a two-term expansion to compute
540 : * extrapolated boundary face values. If this is true, then an implicit system has to be solved. If
541 : * false, then the cell center value will be used as the extrapolated boundary face value
542 : * @param mesh The mesh on which we are computing the gradient
543 : * @return The computed cell gradient
544 : */
545 : template <typename T,
546 : typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
547 : libMesh::VectorValue<T>
548 75 : greenGaussGradient(const FaceArg & face_arg,
549 : const StateArg & state_arg,
550 : const FunctorBase<T> & functor,
551 : const bool two_term_boundary_expansion,
552 : const MooseMesh & mesh)
553 : {
554 : mooseAssert(face_arg.fi, "We should have a face info to compute a face gradient");
555 75 : const auto & fi = *(face_arg.fi);
556 75 : const auto & elem_arg = face_arg.makeElem();
557 75 : const auto & neighbor_arg = face_arg.makeNeighbor();
558 75 : const bool defined_on_elem = functor.hasBlocks(fi.elemSubdomainID());
559 75 : const bool defined_on_neighbor = fi.neighborPtr() && functor.hasBlocks(fi.neighborSubdomainID());
560 :
561 75 : if (defined_on_elem && defined_on_neighbor)
562 : {
563 27 : const auto & value_elem = functor(elem_arg, state_arg);
564 27 : const auto & value_neighbor = functor(neighbor_arg, state_arg);
565 :
566 : // This is the component of the gradient which is parallel to the line connecting
567 : // the cell centers. Therefore, we can use our second order, central difference
568 : // scheme to approximate it.
569 27 : libMesh::VectorValue<T> face_gradient = (value_neighbor - value_elem) / fi.dCNMag() * fi.eCN();
570 :
571 : // We only need nonorthogonal correctors in 2+ dimensions
572 27 : if (mesh.dimension() > 1)
573 : {
574 : // We are using an orthogonal approach for the non-orthogonal correction, for more information
575 : // see Hrvoje Jasak's PhD Thesis (Imperial College, 1996)
576 27 : libMesh::VectorValue<T> interpolated_gradient;
577 :
578 : // Compute the gradients in the two cells on both sides of the face
579 27 : const auto & grad_elem =
580 27 : greenGaussGradient(elem_arg, state_arg, functor, two_term_boundary_expansion, mesh);
581 27 : const auto & grad_neighbor =
582 27 : greenGaussGradient(neighbor_arg, state_arg, functor, two_term_boundary_expansion, mesh);
583 :
584 27 : Moose::FV::interpolate(Moose::FV::InterpMethod::Average,
585 : interpolated_gradient,
586 : grad_elem,
587 : grad_neighbor,
588 : fi,
589 : true);
590 :
591 27 : face_gradient += interpolated_gradient - (interpolated_gradient * fi.eCN()) * fi.eCN();
592 0 : }
593 :
594 27 : return face_gradient;
595 0 : }
596 48 : else if (defined_on_elem)
597 48 : return functor.gradient(elem_arg, state_arg);
598 0 : else if (defined_on_neighbor)
599 0 : return functor.gradient(neighbor_arg, state_arg);
600 : else
601 0 : mooseError("The functor must be defined on one of the sides");
602 : }
603 :
604 : template <typename T>
605 : TensorValue<T>
606 2 : greenGaussGradient(const ElemArg & elem_arg,
607 : const StateArg & state_arg,
608 : const Moose::FunctorBase<libMesh::VectorValue<T>> & functor,
609 : const bool two_term_boundary_expansion,
610 : const MooseMesh & mesh)
611 : {
612 2 : TensorValue<T> ret;
613 8 : for (const auto i : make_range(Moose::dim))
614 : {
615 6 : VectorComponentFunctor<T> scalar_functor(functor, i);
616 0 : const auto row_gradient =
617 6 : greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
618 24 : for (const auto j : make_range(unsigned(Moose::dim)))
619 18 : ret(i, j) = row_gradient(j);
620 : }
621 :
622 2 : return ret;
623 0 : }
624 :
625 : template <typename T>
626 : TensorValue<T>
627 1 : greenGaussGradient(const FaceArg & face_arg,
628 : const StateArg & state_arg,
629 : const Moose::FunctorBase<libMesh::VectorValue<T>> & functor,
630 : const bool two_term_boundary_expansion,
631 : const MooseMesh & mesh)
632 : {
633 1 : TensorValue<T> ret;
634 4 : for (const auto i : make_range(unsigned(Moose::dim)))
635 : {
636 3 : VectorComponentFunctor<T> scalar_functor(functor, i);
637 0 : const auto row_gradient =
638 3 : greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
639 12 : for (const auto j : make_range(unsigned(Moose::dim)))
640 9 : ret(i, j) = row_gradient(j);
641 : }
642 :
643 1 : return ret;
644 0 : }
645 :
646 : template <typename T>
647 : typename Moose::FunctorBase<std::vector<T>>::GradientType
648 80 : greenGaussGradient(const ElemArg & elem_arg,
649 : const StateArg & state_arg,
650 : const Moose::FunctorBase<std::vector<T>> & functor,
651 : const bool two_term_boundary_expansion,
652 : const MooseMesh & mesh)
653 : {
654 : // Determine the size of the container
655 80 : const auto vals = functor(elem_arg, state_arg);
656 : typedef typename Moose::FunctorBase<std::vector<T>>::GradientType GradientType;
657 80 : GradientType ret(vals.size());
658 160 : for (const auto i : index_range(ret))
659 : {
660 : // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
661 : // going to do value type evaluations of the array functor from scalar_functor and we will be
662 : // discarding all the value type evaluations other than the one corresponding to i
663 80 : ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
664 80 : ret[i] =
665 80 : greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
666 : }
667 :
668 160 : return ret;
669 80 : }
670 :
671 : template <typename T>
672 : typename Moose::FunctorBase<std::vector<T>>::GradientType
673 36 : greenGaussGradient(const FaceArg & face_arg,
674 : const StateArg & state_arg,
675 : const Moose::FunctorBase<std::vector<T>> & functor,
676 : const bool two_term_boundary_expansion,
677 : const MooseMesh & mesh)
678 : {
679 : // Determine the size of the container
680 36 : const auto vals = functor(face_arg, state_arg);
681 : typedef typename Moose::FunctorBase<std::vector<T>>::GradientType GradientType;
682 36 : GradientType ret(vals.size());
683 72 : for (const auto i : index_range(ret))
684 : {
685 : // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
686 : // going to do value type evaluations of the array functor from scalar_functor and we will be
687 : // discarding all the value type evaluations other than the one corresponding to i
688 36 : ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
689 36 : ret[i] =
690 36 : greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
691 : }
692 :
693 72 : return ret;
694 36 : }
695 :
696 : template <typename T, std::size_t N>
697 : typename Moose::FunctorBase<std::array<T, N>>::GradientType
698 76 : greenGaussGradient(const ElemArg & elem_arg,
699 : const StateArg & state_arg,
700 : const Moose::FunctorBase<std::array<T, N>> & functor,
701 : const bool two_term_boundary_expansion,
702 : const MooseMesh & mesh)
703 : {
704 : typedef typename Moose::FunctorBase<std::array<T, N>>::GradientType GradientType;
705 76 : GradientType ret;
706 152 : for (const auto i : make_range(N))
707 : {
708 : // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
709 : // going to do value type evaluations of the array functor from scalar_functor and we will be
710 : // discarding all the value type evaluations other than the one corresponding to i
711 76 : ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
712 76 : ret[i] =
713 76 : greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
714 : }
715 :
716 76 : return ret;
717 : }
718 :
719 : template <typename T, std::size_t N>
720 : typename Moose::FunctorBase<std::array<T, N>>::GradientType
721 36 : greenGaussGradient(const FaceArg & face_arg,
722 : const StateArg & state_arg,
723 : const Moose::FunctorBase<std::array<T, N>> & functor,
724 : const bool two_term_boundary_expansion,
725 : const MooseMesh & mesh)
726 : {
727 : typedef typename Moose::FunctorBase<std::array<T, N>>::GradientType GradientType;
728 36 : GradientType ret;
729 72 : for (const auto i : make_range(N))
730 : {
731 : // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
732 : // going to do value type evaluations of the array functor from scalar_functor and we will be
733 : // discarding all the value type evaluations other than the one corresponding to i
734 36 : ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
735 36 : ret[i] =
736 36 : greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
737 : }
738 :
739 36 : return ret;
740 : }
741 : }
742 : }
|
