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 27200 : return libMesh::outer_product(a, b);
40 : }
41 : }
42 :
43 : template <typename T, typename Map>
44 : typename FaceCenteredMapFunctor<T, Map>::ValueType
45 6800 : 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 6800 : 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 6800 : DenseMatrix<PrimitiveType> n_x_Sf(dim, dim);
79 6800 : DenseVector<PrimitiveType> sum_normal_flux(dim);
80 :
81 6800 : const Elem * const elem = elem_arg.elem;
82 :
83 34000 : 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 27200 : const bool elem_has_fi = Moose::FV::elemHasFaceInfo(*elem, neighbor);
90 40200 : const FaceInfo * const fi = _mesh.faceInfo(
91 13000 : elem_has_fi ? elem : neighbor, elem_has_fi ? side : neighbor->which_neighbor_am_i(elem));
92 27200 : const Point & normal = elem_has_fi ? fi->normal() : Point(-fi->normal());
93 :
94 : const Point area_vector = normal * fi->faceArea();
95 27200 : 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 81600 : for (const auto i : make_range(dim))
101 : {
102 54400 : sum_normal_flux(i) += flux_contrib(i);
103 163200 : for (const auto j : make_range(dim))
104 108800 : 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 6800 : DenseVector<PrimitiveType> dense_result(dim);
111 6800 : n_x_Sf.cholesky_solve(sum_normal_flux, dense_result);
112 :
113 : ValueType result;
114 20400 : for (const auto i : make_range(dim))
115 13600 : result(i) = dense_result(i);
116 :
117 13600 : return result;
118 6800 : }
119 : else
120 0 : mooseError("Cell center reconstruction is not implemented!");
121 : }
122 :
123 : template <typename T, typename Map>
124 : typename FaceCenteredMapFunctor<T, Map>::ValueType
125 135621260 : FaceCenteredMapFunctor<T, Map>::evaluate(const FaceArg & face, const StateArg &) const
126 : {
127 135621260 : return this->evaluate(face.fi);
128 : }
129 :
130 : template <typename T, typename Map>
131 : typename FaceCenteredMapFunctor<T, Map>::ValueType
132 203461643 : FaceCenteredMapFunctor<T, Map>::evaluate(const FaceInfo * fi) const
133 : {
134 : try
135 : {
136 203461643 : return libmesh_map_find(*this, fi->id());
137 : }
138 2 : catch (libMesh::LogicError &)
139 : {
140 2 : if (!_sub_ids.empty() && !_sub_ids.count(fi->elem().subdomain_id()))
141 : {
142 1 : if (fi->neighborPtr() && !_sub_ids.count(fi->neighborPtr()->subdomain_id()))
143 1 : mooseError("Attempted to evaluate FaceCenteredMapFunctor '",
144 : this->functorName(),
145 : "' with an element subdomain id of '",
146 1 : fi->elem().subdomain_id(),
147 5 : fi->neighborPtr() ? " or neighbor subdomain id of '" +
148 : std::to_string(fi->neighborPtr()->subdomain_id()) + "'"
149 : : "",
150 : "' but that subdomain id is not one of the subdomain ids the functor is "
151 : "restricted to.");
152 : }
153 : else
154 1 : mooseError("Attempted access into FaceCenteredMapFunctor '",
155 : this->functorName(),
156 : "' with a key that does not yet exist in the map. Make sure to fill your "
157 : "FaceCenteredMapFunctor for all elements you will attempt to access later.");
158 :
159 : return typename FaceCenteredMapFunctor<T, Map>::ValueType();
160 : }
161 : }
162 :
163 : template class FaceCenteredMapFunctor<ADRealVectorValue,
164 : std::unordered_map<dof_id_type, ADRealVectorValue>>;
165 : template class FaceCenteredMapFunctor<RealVectorValue,
166 : std::unordered_map<dof_id_type, RealVectorValue>>;
167 : template class FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>;
|
