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 "FaceCenteredMapFunctor.h"
11 : #include "FaceInfo.h"
12 : #include "FVUtils.h"
13 : #include "libmesh/compare_types.h"
14 : #include "libmesh/type_tensor.h"
15 : #include "libmesh/tensor_tools.h"
16 : #include "libmesh/dense_matrix.h"
17 : #include "libmesh/elem.h"
18 : #include "libmesh/point.h"
19 :
20 : using namespace libMesh;
21 :
22 : namespace Moose
23 : {
24 : template <typename T, typename T2, typename std::enable_if<ScalarTraits<T>::value, int>::type = 0>
25 : inline TypeVector<typename CompareTypes<T, T2>::supertype>
26 : outer_product(const T & a, const TypeVector<T2> & b)
27 : {
28 : TypeVector<typename CompareTypes<T, T2>::supertype> ret;
29 : for (const auto i : make_range(Moose::dim))
30 : ret(i) = a * b(i);
31 :
32 : return ret;
33 : }
34 :
35 : template <typename T, typename T2>
36 : inline TypeTensor<typename CompareTypes<T, T2>::supertype>
37 : outer_product(const TypeVector<T> & a, const TypeVector<T2> & b)
38 : {
39 16320 : return libMesh::outer_product(a, b);
40 : }
41 : }
42 :
43 : template <typename T, typename Map>
44 : typename FaceCenteredMapFunctor<T, Map>::ValueType
45 4080 : FaceCenteredMapFunctor<T, Map>::evaluate(const ElemArg & elem_arg, const StateArg &) const
46 : {
47 : // The following reconstruction is based on Weller's method. For more information on this,
48 : // we recommend:
49 : //
50 : // Weller, Hilary. "Non-orthogonal version of the arbitrary polygonal C-grid and a new diamond
51 : // grid." Geoscientific Model Development 7.3 (2014): 779-797.
52 : //
53 : // and
54 : //
55 : // Aguerre, Horacio J., et al. "An oscillation-free flow solver based on flux reconstruction."
56 : // Journal of Computational Physics 365 (2018): 135-148.
57 : //
58 : // This basically reconstructs the cell value based on flux values as follows:
59 : //
60 : // $\left( \sum_f n_f \outer S_f \right)^{-1} \sum_f (\phi_f \cdot S_f)n_f$
61 : //
62 : // where $S_f$ is the surface normal vector, $n_f$ is the unit surface vector and $\phi_f$ is a
63 : // vector value on the field. Hence the restriction to vector values.
64 :
65 : using ValueType = typename FaceCenteredMapFunctor<T, Map>::ValueType;
66 :
67 : // Primitive type one rank below the type of the stored data (T). If we are storing a rank two
68 : // tensor, this is a vector, if we are storing a vector this is just a number
69 : using PrimitiveType = typename MetaPhysicL::ReplaceAlgebraicType<
70 : T,
71 : typename TensorTools::DecrementRank<typename MetaPhysicL::ValueType<T>::type>::type>::type;
72 :
73 : if constexpr (libMesh::TensorTools::TensorTraits<ValueType>::rank == 1)
74 : {
75 4080 : const auto dim = _mesh.dimension();
76 : // The reason why these are DenseVector/Matrix is that when the mesh dimension is lower than 3,
77 : // we get singular matrixes if we try to invert TensorValues.
78 4080 : DenseMatrix<PrimitiveType> n_x_Sf(dim, dim);
79 4080 : DenseVector<PrimitiveType> sum_normal_flux(dim);
80 :
81 4080 : const Elem * const elem = elem_arg.elem;
82 :
83 20400 : for (const auto side : make_range(elem->n_sides()))
84 : {
85 : const Elem * const neighbor = elem->neighbor_ptr(side);
86 :
87 : // We need to check if the faceinfo belongs to the element or the neighbor. Based on that we
88 : // query the faceinfo and adjust the normal to point outward of the current cell
89 16320 : const bool elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
90 24120 : const FaceInfo * const fi = _mesh.faceInfo(
91 7800 : elem_has_fi ? elem : neighbor, elem_has_fi ? side : neighbor->which_neighbor_am_i(elem));
92 16320 : const Point & normal = elem_has_fi ? fi->normal() : Point(-fi->normal());
93 :
94 : const Point area_vector = normal * fi->faceArea();
95 16320 : const ValueType face_value = this->evaluate(fi);
96 :
97 : const auto product = Moose::outer_product(normal, area_vector);
98 :
99 0 : const auto flux_contrib = normal * (face_value * area_vector);
100 48960 : for (const auto i : make_range(dim))
101 : {
102 32640 : sum_normal_flux(i) += flux_contrib(i);
103 97920 : for (const auto j : make_range(dim))
104 65280 : n_x_Sf(i, j) += product(i, j);
105 : }
106 : }
107 :
108 : // We do the inversion of the surface vector matrix here. It is symmetric
109 : // and small so we can do it using a Cholesky decomposition.
110 4080 : DenseVector<PrimitiveType> dense_result(dim);
111 4080 : n_x_Sf.cholesky_solve(sum_normal_flux, dense_result);
112 :
113 : ValueType result;
114 12240 : for (const auto i : make_range(dim))
115 8160 : result(i) = dense_result(i);
116 :
117 8160 : return result;
118 4080 : }
119 : else
120 : {
121 : (void)elem_arg; // WE do this because the GCC min complains that it is not used and
122 : // [[maybe_unused]] doesnt work for GCC min
123 0 : mooseError("Cell center reconstruction is not implemented!");
124 : }
125 : }
126 :
127 : template <typename T, typename Map>
128 : typename FaceCenteredMapFunctor<T, Map>::ValueType
129 116314204 : FaceCenteredMapFunctor<T, Map>::evaluate(const FaceArg & face, const StateArg &) const
130 : {
131 116314204 : return this->evaluate(face.fi);
132 : }
133 :
134 : template <typename T, typename Map>
135 : typename FaceCenteredMapFunctor<T, Map>::ValueType
136 183520557 : FaceCenteredMapFunctor<T, Map>::evaluate(const FaceInfo * fi) const
137 : {
138 183520557 : auto it = this->find(fi->id());
139 183520557 : if (it == this->end())
140 : {
141 2 : if (!_sub_ids.empty() && !_sub_ids.count(fi->elem().subdomain_id()))
142 : {
143 1 : if (fi->neighborPtr() && !_sub_ids.count(fi->neighborPtr()->subdomain_id()))
144 1 : mooseError("Attempted to evaluate FaceCenteredMapFunctor '",
145 : this->functorName(),
146 : "' with an element subdomain id of '",
147 1 : fi->elem().subdomain_id(),
148 3 : fi->neighborPtr() ? " or neighbor subdomain id of '" +
149 : std::to_string(fi->neighborPtr()->subdomain_id()) + "'"
150 : : "",
151 : "' but that subdomain id is not one of the subdomain ids the functor is "
152 : "restricted to.");
153 : }
154 : else
155 1 : mooseError("Attempted access into FaceCenteredMapFunctor '",
156 : this->functorName(),
157 : "' with a key that does not yet exist in the map. Make sure to fill your "
158 : "FaceCenteredMapFunctor for all elements you will attempt to access later.");
159 :
160 0 : return typename FaceCenteredMapFunctor<T, Map>::ValueType();
161 : }
162 :
163 183520555 : return it->second;
164 : }
165 :
166 : template class FaceCenteredMapFunctor<ADRealVectorValue,
167 : std::unordered_map<dof_id_type, ADRealVectorValue>>;
168 : template class FaceCenteredMapFunctor<RealVectorValue,
169 : std::unordered_map<dof_id_type, RealVectorValue>>;
170 : template class FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>;
|
