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 "FunctorSmoother.h"
11 : #include "MooseUtils.h"
12 : #include "RelationshipManager.h"
13 :
14 : #include "metaphysicl/raw_type.h"
15 :
16 : registerMooseObject("MooseApp", FunctorSmoother);
17 :
18 : template <typename T>
19 : InputParameters
20 14346 : FunctorSmootherTempl<T>::validParams()
21 : {
22 14346 : InputParameters params = FunctorMaterial::validParams();
23 14346 : params.addClassDescription("Creates smoother functor(s) using various averaging techniques");
24 14346 : params.addParam<std::vector<MooseFunctorName>>("functors_in",
25 : "The name(s) of the functors to smooth");
26 14346 : params.addParam<std::vector<MooseFunctorName>>("functors_out",
27 : "The name(s) of the smooth output functors");
28 14346 : MooseEnum smoothing_techniques("face_average layered_elem_average remove_checkerboard");
29 14346 : params.addParam<MooseEnum>(
30 : "smoothing_technique", smoothing_techniques, "How to smooth the functor");
31 :
32 28692 : return params;
33 14346 : }
34 :
35 : template <typename T>
36 42 : FunctorSmootherTempl<T>::FunctorSmootherTempl(const InputParameters & parameters)
37 : : FunctorMaterial(parameters),
38 42 : _functors_in(getParam<std::vector<MooseFunctorName>>("functors_in")),
39 42 : _functors_out(getParam<std::vector<MooseFunctorName>>("functors_out")),
40 84 : _smoothing_technique(getParam<MooseEnum>("smoothing_technique"))
41 : {
42 42 : if (_tid == 0)
43 : {
44 : // We need one layer of ghosting at least to get neighbor values
45 39 : auto & factory = _app.getFactory();
46 39 : auto rm_params = factory.getValidParams("ElementSideNeighborLayers");
47 39 : rm_params.template set<std::string>("for_whom") = name();
48 39 : rm_params.template set<MooseMesh *>("mesh") = &const_cast<MooseMesh &>(_mesh);
49 39 : rm_params.template set<Moose::RelationshipManagerType>("rm_type") =
50 39 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
51 : Moose::RelationshipManagerType::COUPLING;
52 39 : rm_params.template set<unsigned short>("layers") = 1;
53 39 : rm_params.template set<bool>("use_point_neighbors") = false;
54 39 : rm_params.template set<bool>("use_displaced_mesh") =
55 39 : parameters.template get<bool>("use_displaced_mesh");
56 : mooseAssert(rm_params.areAllRequiredParamsValid(),
57 : "All relationship manager parameters should be valid.");
58 39 : auto rm_obj = factory.template create<RelationshipManager>(
59 39 : "ElementSideNeighborLayers", name() + "_functor_smoothing", rm_params);
60 :
61 : // Delete the resources created on behalf of the RM if it ends up not being added to the
62 : // App.
63 39 : if (!_app.addRelationshipManager(rm_obj))
64 39 : factory.releaseSharedObjects(*rm_obj);
65 39 : }
66 :
67 42 : const std::set<ExecFlagType> clearance_schedule(_execute_enum.begin(), _execute_enum.end());
68 :
69 84 : for (const auto i : index_range(_functors_in))
70 : {
71 : // This potentially upcasts functors from non-AD to AD
72 42 : const auto & functor_in = getFunctor<typename Moose::ADType<T>::type>(_functors_in[i]);
73 :
74 : // We always add the AD type
75 42 : addFunctorProperty<typename Moose::ADType<T>::type>(
76 : _functors_out[i],
77 19737 : [this, &functor_in](const auto & r, const auto & t) -> typename Moose::ADType<T>::type
78 : {
79 2700 : typename Moose::ADType<T>::type base_value = functor_in(r, t);
80 2700 : typename Moose::ADType<T>::type average = 0;
81 :
82 2700 : const Elem * r_elem = nullptr;
83 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
84 2700 : r_elem = r.elem;
85 :
86 : // Handle the single isolated element case
87 2700 : if (r_elem && r_elem->n_neighbors() == 0)
88 0 : return base_value;
89 :
90 : // If these techniques are of interest elsewhere, they should be moved to a
91 : // FunctorAveragingUtils file and namespaced
92 2700 : if (_smoothing_technique == FACE_AVERAGE)
93 : {
94 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
95 : {
96 900 : unsigned int n_faces = 0;
97 4500 : for (const auto side_index : r_elem->side_index_range())
98 : {
99 3600 : auto fi = _mesh.faceInfo(r_elem, side_index);
100 3600 : if (!fi)
101 : fi =
102 891 : _mesh.faceInfo(r_elem->neighbor_ptr(side_index),
103 : r_elem->neighbor_ptr(side_index)->which_neighbor_am_i(r_elem));
104 3600 : Moose::FaceArg face_arg{
105 : fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
106 3600 : if (face_arg.fi)
107 : {
108 3600 : average += functor_in(face_arg, t);
109 3600 : n_faces++;
110 : }
111 : }
112 900 : average /= n_faces;
113 : }
114 : else
115 : // a conceptually simple option here would be to use the smoothed version of the
116 : // ElemArg to compute all the other functor arguments. Maybe with the same smoothing
117 : // on the gradient calls
118 0 : mooseError("Face averaging smoothing has only been defined for the ElemArg functor "
119 : "argument, not for ",
120 : MooseUtils::prettyCppType(&r),
121 : ". Please contact a MOOSE developer or implement it yourself.");
122 : }
123 1800 : else if (_smoothing_technique == LAYERED_AVERAGE)
124 : {
125 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
126 : {
127 900 : unsigned int n_neighbors = 0;
128 4500 : for (const auto neighbor : r_elem->neighbor_ptr_range())
129 3600 : if (neighbor)
130 : {
131 1782 : n_neighbors++;
132 1782 : average += functor_in(Moose::ElemArg{neighbor, false}, t);
133 : }
134 900 : average /= n_neighbors;
135 : }
136 : else
137 0 : mooseError("Element layered averaging smoothing has only been defined for the "
138 : "ElemArg functor argument, not for ",
139 : MooseUtils::prettyCppType(&r),
140 : ". Please contact a MOOSE developer or implement it yourself.");
141 : }
142 : else
143 : {
144 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
145 : {
146 : // We average the local value with the average of the two values furthest away
147 : // among the neighbors. For a checkerboarding but overall linear evolution, this will
148 : // compute the average of the "high" and "low" profiles
149 900 : average += base_value / 2;
150 :
151 : // Find the neighbor with the furthest value
152 900 : typename Moose::ADType<T>::type extreme_value = 0;
153 900 : typename Moose::ADType<T>::type delta = 0;
154 :
155 900 : unsigned int furthest_one = 0;
156 900 : unsigned int num_neighbors = 0;
157 4500 : for (const auto side_index : r_elem->side_index_range())
158 : {
159 3600 : auto neighbor = r_elem->neighbor_ptr(side_index);
160 3600 : if (neighbor)
161 : {
162 1782 : num_neighbors++;
163 1782 : auto neighbor_value = functor_in(Moose::ElemArg{neighbor, false}, t);
164 1782 : if (std::abs(neighbor_value - base_value) > delta)
165 : {
166 900 : furthest_one = side_index;
167 900 : extreme_value = neighbor_value;
168 900 : delta = std::abs(neighbor_value - base_value);
169 : }
170 1782 : }
171 : }
172 :
173 : // We're on a boundary in 1D, or maybe an odd shape corner in 2D
174 900 : if (num_neighbors == 1)
175 : {
176 18 : average += extreme_value / 2;
177 18 : return average;
178 : }
179 : else
180 882 : average += extreme_value / 4;
181 :
182 : // Get the value from the neighbor opposite the furthest-value neighbor
183 : try
184 : {
185 882 : auto opposite_side = r_elem->opposite_side(furthest_one);
186 882 : auto neighbor = r_elem->neighbor_ptr(opposite_side);
187 882 : if (neighbor)
188 : {
189 882 : average += functor_in(Moose::ElemArg{neighbor, false}, t) / 4;
190 : }
191 : else
192 : // We're probably at a boundary
193 0 : average += extreme_value / 4;
194 : }
195 : // if there is no opposite side (for example a triangle)
196 0 : catch (libMesh::LogicError & e)
197 : {
198 : // find the second furthest
199 0 : delta = 0;
200 0 : typename Moose::ADType<T>::type second_extreme_value = 0;
201 :
202 0 : for (const auto side_index : r_elem->side_index_range())
203 : {
204 0 : auto neighbor = r_elem->neighbor_ptr(side_index);
205 0 : auto neighbor_value =
206 0 : neighbor ? functor_in(Moose::ElemArg{neighbor, false}, t) : base_value;
207 0 : if (std::abs(neighbor_value - base_value) > delta && side_index != furthest_one)
208 : {
209 0 : second_extreme_value = neighbor_value;
210 0 : delta = std::abs(neighbor_value - base_value);
211 : }
212 : }
213 0 : average += second_extreme_value / 4;
214 0 : }
215 918 : }
216 : else
217 0 : mooseError("Checkerboard removal smoothing has only been defined for the "
218 : "ElemArg functor argument, not for ",
219 : MooseUtils::prettyCppType(&r),
220 : ". Please contact a MOOSE developer or implement it yourself.");
221 : }
222 2682 : return average;
223 2700 : },
224 : clearance_schedule);
225 : }
226 42 : }
227 :
228 : template class FunctorSmootherTempl<Real>;
|
