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