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 14340 : FunctorSmootherTempl<T>::validParams()
21 : {
22 14340 : InputParameters params = FunctorMaterial::validParams();
23 14340 : params.addClassDescription("Creates smoother functor(s) using various averaging techniques");
24 14340 : params.addParam<std::vector<MooseFunctorName>>("functors_in",
25 : "The name(s) of the functors to smooth");
26 14340 : params.addParam<std::vector<MooseFunctorName>>("functors_out",
27 : "The name(s) of the smooth output functors");
28 14340 : MooseEnum smoothing_techniques("face_average layered_elem_average remove_checkerboard");
29 14340 : params.addParam<MooseEnum>(
30 : "smoothing_technique", smoothing_techniques, "How to smooth the functor");
31 :
32 28680 : return params;
33 14340 : }
34 :
35 : template <typename T>
36 39 : FunctorSmootherTempl<T>::FunctorSmootherTempl(const InputParameters & parameters)
37 : : FunctorMaterial(parameters),
38 39 : _functors_in(getParam<std::vector<MooseFunctorName>>("functors_in")),
39 39 : _functors_out(getParam<std::vector<MooseFunctorName>>("functors_out")),
40 78 : _smoothing_technique(getParam<MooseEnum>("smoothing_technique"))
41 : {
42 39 : if (_tid == 0)
43 : {
44 : // We need one layer of ghosting at least to get neighbor values
45 36 : auto & factory = _app.getFactory();
46 36 : auto rm_params = factory.getValidParams("ElementSideNeighborLayers");
47 36 : rm_params.template set<std::string>("for_whom") = name();
48 36 : rm_params.template set<MooseMesh *>("mesh") = &const_cast<MooseMesh &>(_mesh);
49 36 : rm_params.template set<Moose::RelationshipManagerType>("rm_type") =
50 36 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
51 : Moose::RelationshipManagerType::COUPLING;
52 36 : rm_params.template set<unsigned short>("layers") = 1;
53 36 : rm_params.template set<bool>("use_point_neighbors") = false;
54 36 : rm_params.template set<bool>("use_displaced_mesh") =
55 36 : parameters.template get<bool>("use_displaced_mesh");
56 : mooseAssert(rm_params.areAllRequiredParamsValid(),
57 : "All relationship manager parameters should be valid.");
58 36 : auto rm_obj = factory.template create<RelationshipManager>(
59 36 : "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 36 : if (!_app.addRelationshipManager(rm_obj))
64 36 : factory.releaseSharedObjects(*rm_obj);
65 36 : }
66 :
67 39 : const std::set<ExecFlagType> clearance_schedule(_execute_enum.begin(), _execute_enum.end());
68 :
69 78 : for (const auto i : index_range(_functors_in))
70 : {
71 : // This potentially upcasts functors from non-AD to AD
72 39 : const auto & functor_in = getFunctor<typename Moose::ADType<T>::type>(_functors_in[i]);
73 :
74 : // We always add the AD type
75 39 : addFunctorProperty<typename Moose::ADType<T>::type>(
76 : _functors_out[i],
77 17544 : [this, &functor_in](const auto & r, const auto & t) -> typename Moose::ADType<T>::type
78 : {
79 2400 : typename Moose::ADType<T>::type base_value = functor_in(r, t);
80 2400 : typename Moose::ADType<T>::type average = 0;
81 :
82 2400 : const Elem * r_elem = nullptr;
83 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
84 2400 : r_elem = r.elem;
85 :
86 : // Handle the single isolated element case
87 2400 : 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 2400 : if (_smoothing_technique == FACE_AVERAGE)
93 : {
94 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
95 : {
96 800 : unsigned int n_faces = 0;
97 4000 : for (const auto side_index : r_elem->side_index_range())
98 : {
99 3200 : auto fi = _mesh.faceInfo(r_elem, side_index);
100 3200 : if (!fi)
101 : fi =
102 792 : _mesh.faceInfo(r_elem->neighbor_ptr(side_index),
103 : r_elem->neighbor_ptr(side_index)->which_neighbor_am_i(r_elem));
104 3200 : Moose::FaceArg face_arg{
105 : fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
106 3200 : if (face_arg.fi)
107 : {
108 3200 : average += functor_in(face_arg, t);
109 3200 : n_faces++;
110 : }
111 : }
112 800 : 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 1600 : else if (_smoothing_technique == LAYERED_AVERAGE)
124 : {
125 : if constexpr (std::is_same_v<const Moose::ElemArg &, decltype(r)>)
126 : {
127 800 : unsigned int n_neighbors = 0;
128 4000 : for (const auto neighbor : r_elem->neighbor_ptr_range())
129 3200 : if (neighbor)
130 : {
131 1584 : n_neighbors++;
132 1584 : average += functor_in(Moose::ElemArg{neighbor, false}, t);
133 : }
134 800 : 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 800 : average += base_value / 2;
150 :
151 : // Find the neighbor with the furthest value
152 800 : typename Moose::ADType<T>::type extreme_value = 0;
153 800 : typename Moose::ADType<T>::type delta = 0;
154 :
155 800 : unsigned int furthest_one = 0;
156 800 : unsigned int num_neighbors = 0;
157 4000 : for (const auto side_index : r_elem->side_index_range())
158 : {
159 3200 : auto neighbor = r_elem->neighbor_ptr(side_index);
160 3200 : if (neighbor)
161 : {
162 1584 : num_neighbors++;
163 1584 : auto neighbor_value = functor_in(Moose::ElemArg{neighbor, false}, t);
164 1584 : if (std::abs(neighbor_value - base_value) > delta)
165 : {
166 800 : furthest_one = side_index;
167 800 : extreme_value = neighbor_value;
168 800 : delta = std::abs(neighbor_value - base_value);
169 : }
170 1584 : }
171 : }
172 :
173 : // We're on a boundary in 1D, or maybe an odd shape corner in 2D
174 800 : if (num_neighbors == 1)
175 : {
176 16 : average += extreme_value / 2;
177 16 : return average;
178 : }
179 : else
180 784 : average += extreme_value / 4;
181 :
182 : // Get the value from the neighbor opposite the furthest-value neighbor
183 : try
184 : {
185 784 : auto opposite_side = r_elem->opposite_side(furthest_one);
186 784 : auto neighbor = r_elem->neighbor_ptr(opposite_side);
187 784 : if (neighbor)
188 : {
189 784 : 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 816 : }
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 2384 : return average;
223 2400 : },
224 : clearance_schedule);
225 : }
226 39 : }
227 :
228 : template class FunctorSmootherTempl<Real>;
|
