https://mooseframework.inl.gov
Functions
TestReconstruction.C File Reference

Go to the source code of this file.

Functions

void testReconstruction (const Moose::CoordinateSystemType coord_type)
 
 TEST (TestReconstruction, Cartesian)
 
 TEST (TestReconstruction, Cylindrical)
 

Function Documentation

◆ TEST() [1/2]

TEST ( TestReconstruction  ,
Cartesian   
)

Definition at line 242 of file TestReconstruction.C.

void testReconstruction(const Moose::CoordinateSystemType coord_type)

◆ TEST() [2/2]

TEST ( TestReconstruction  ,
Cylindrical   
)

Definition at line 244 of file TestReconstruction.C.

void testReconstruction(const Moose::CoordinateSystemType coord_type)

◆ testReconstruction()

void testReconstruction ( const Moose::CoordinateSystemType  coord_type)

Definition at line 41 of file TestReconstruction.C.

Referenced by TEST().

42 {
43  MultiMooseEnum coord_type_enum("XYZ RZ RSPHERICAL", "XYZ");
44  coord_type_enum = (coord_type == Moose::COORD_XYZ) ? "XYZ" : "RZ";
45 
46  std::vector<unsigned int> num_elem = {64, 128, 256};
47  std::vector<Real> weller_errors;
48  std::vector<Real> moukalled_errors;
49  std::vector<Real> tano_errors;
50  std::vector<Real> tano_errors_twice;
51  std::vector<Real> h(num_elem.size());
52  for (const auto i : index_range(num_elem))
53  h[i] = 1. / num_elem[i];
54 
55  for (const auto i : index_range(num_elem))
56  {
57  const auto nx = num_elem[i];
58  auto app = AppFactory::create("NavierStokesUnitApp");
59  auto * factory = &app->getFactory();
60  std::string mesh_type = "MeshGeneratorMesh";
61 
62  std::shared_ptr<MeshGeneratorMesh> mesh;
63  {
64  InputParameters params = factory->getValidParams(mesh_type);
65  mesh = factory->create<MeshGeneratorMesh>(mesh_type, "moose_mesh", params);
66  }
67 
68  app->actionWarehouse().mesh() = mesh;
69 
70  {
71  std::unique_ptr<MeshBase> lm_mesh;
72  InputParameters params = factory->getValidParams("GeneratedMeshGenerator");
73  params.set<unsigned int>("nx") = nx;
74  params.set<unsigned int>("ny") = nx;
75  params.set<MooseEnum>("dim") = "2";
76  auto mesh_gen =
77  factory->create<GeneratedMeshGenerator>("GeneratedMeshGenerator", "mesh_gen", params);
78  lm_mesh = mesh_gen->generate();
79  mesh->setMeshBase(std::move(lm_mesh));
80  }
81 
82  mesh->prepare(nullptr);
83  mesh->setCoordSystem({}, coord_type_enum);
84  mooseAssert(mesh->getAxisymmetricRadialCoord() == 0,
85  "This should be 0 because we haven't set anything.");
86  const auto & all_fi = mesh->allFaceInfo();
87  mesh->buildFiniteVolumeInfo();
88  mesh->computeFiniteVolumeCoords();
89  std::vector<const FaceInfo *> faces(all_fi.size());
90  for (const auto i : index_range(all_fi))
91  faces[i] = &all_fi[i];
92 
93  auto & lm_mesh = mesh->getMesh();
94 
96  *mesh, "u", /*extrapoalted_bd=*/true);
97  for (auto * const elem : lm_mesh.active_element_ptr_range())
98  {
99  const auto centroid = elem->vertex_average();
100  const auto value = RealVectorValue(-std::sin(centroid(0)) * std::cos(centroid(1)),
101  std::cos(centroid(0)) * std::sin(centroid(1)));
102  u[elem->id()] = value;
103  }
104 
106  up_weller(*mesh, "up_weller", /*extrapoalted_bd=*/true);
108  up_moukalled(*mesh, "up_moukalled", /*extrapoalted_bd=*/true);
110  up_tano(*mesh, "up_tano", /*extrapoalted_bd=*/true);
111 
112  for (const auto & fi : all_fi)
113  {
114  auto moukalled_reconstruct = [&fi](auto & functor, auto & container)
115  {
116  auto face = Moose::FaceArg(
117  {&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr});
118  const RealVectorValue uf(functor(face, Moose::currentState()));
119  const Point surface_vector = fi.normal() * fi.faceArea();
120  auto product = (uf * fi.dCN()) * surface_vector;
121 
122  container[fi.elem().id()] += product * fi.gC() / fi.elemVolume();
123  if (fi.neighborPtr())
124  container[fi.neighbor().id()] +=
125  std::move(product) * (1. - fi.gC()) / fi.neighborVolume();
126  };
127 
128  moukalled_reconstruct(u, up_moukalled);
129  }
130 
131  Moose::FV::interpolateReconstruct(up_weller, u, 1, true, faces, Moose::currentState());
132  Moose::FV::interpolateReconstruct(up_tano, u, 1, false, faces, Moose::currentState());
133 
134  Real error = 0;
135  Real weller_error = 0;
136  Real tano_error = 0;
137  Real moukalled_error = 0;
138  const auto current_h = h[i];
139  for (auto * const elem : lm_mesh.active_element_ptr_range())
140  {
141  const auto elem_id = elem->id();
142  auto elem_arg = Moose::ElemArg{elem, false};
143  const RealVectorValue analytic(u(elem_arg, Moose::currentState()));
144 
145  auto compute_elem_error =
146  [elem_id, current_h, elem, &analytic](auto & container, auto & error)
147  {
148  auto & current = libmesh_map_find(container, elem_id);
149  const auto diff = analytic - current;
150  error += diff * diff * current_h * current_h;
151 
152  // Test CellCenteredMapFunctor ElemPointArg overload
153  const auto elem_point_eval = container(
154  Moose::ElemPointArg({elem, elem->vertex_average(), false}), Moose::currentState());
155  for (const auto d : make_range(Moose::dim))
156  EXPECT_TRUE(MooseUtils::absoluteFuzzyEqual(current(d), elem_point_eval(d)));
157  };
158 
159  compute_elem_error(up_weller, weller_error);
160  compute_elem_error(up_moukalled, moukalled_error);
161  compute_elem_error(up_tano, tano_error);
162  }
163  error = std::sqrt(error);
164  weller_error = std::sqrt(weller_error);
165  moukalled_error = std::sqrt(moukalled_error);
166  tano_error = std::sqrt(tano_error);
167  weller_errors.push_back(weller_error);
168  moukalled_errors.push_back(moukalled_error);
169  tano_errors.push_back(tano_error);
170 
171  up_tano.clear();
172  Moose::FV::interpolateReconstruct(up_tano, u, 2, false, faces, Moose::currentState());
173 
174  tano_error = 0;
175  for (auto * const elem : lm_mesh.active_element_ptr_range())
176  {
177  const auto elem_id = elem->id();
178  const auto elem_arg = Moose::ElemArg{elem, false};
179  const RealVectorValue analytic(u(elem_arg, Moose::currentState()));
180 
181  auto compute_elem_error = [elem_id, current_h, &analytic](auto & container, auto & error)
182  {
183  auto & current = libmesh_map_find(container, elem_id);
184  const auto diff = analytic - current;
185  error += diff * diff * current_h * current_h;
186  };
187  compute_elem_error(up_tano, tano_error);
188  }
189  tano_error = std::sqrt(tano_error);
190  tano_errors_twice.push_back(tano_error);
191 
193  *mesh, "not_restricted", /*extrapoalted_bd=*/true);
194  try
195  {
196  unrestricted_error_test(Moose::ElemArg({lm_mesh.elem_ptr(0), false}), Moose::currentState());
197  EXPECT_TRUE(false);
198  }
199  catch (std::runtime_error & e)
200  {
201  EXPECT_TRUE(std::string(e.what()).find("not_restricted") != std::string::npos);
202  EXPECT_TRUE(std::string(e.what()).find("Make sure to fill") != std::string::npos);
203  }
204 
206  *mesh, {1}, "is_restricted", /*extrapoalted_bd=*/true);
207  try
208  {
209  restricted_error_test(Moose::ElemArg({lm_mesh.elem_ptr(0), false}), Moose::currentState());
210  EXPECT_TRUE(false);
211  }
212  catch (std::runtime_error & e)
213  {
214  EXPECT_TRUE(std::string(e.what()).find("is_restricted") != std::string::npos);
215  EXPECT_TRUE(std::string(e.what()).find("0") != std::string::npos);
216  EXPECT_TRUE(
217  std::string(e.what()).find(
218  "that subdomain id is not one of the subdomain ids the functor is restricted to") !=
219  std::string::npos);
220  }
221  }
222 
223  std::for_each(h.begin(), h.end(), [](Real & h_elem) { h_elem = std::log(h_elem); });
224 
225  auto expect_errors = [&h](auto & errors_arg, Real expected_error)
226  {
227  std::for_each(
228  errors_arg.begin(), errors_arg.end(), [](Real & error) { error = std::log(error); });
229  PolynomialFit fit(h, errors_arg, 1);
230  fit.generate();
231 
232  const auto & coeffs = fit.getCoefficients();
233  EXPECT_NEAR(coeffs[1], expected_error, .05);
234  };
235 
236  expect_errors(weller_errors, coord_type == Moose::COORD_RZ ? 1.5 : 2);
237  expect_errors(moukalled_errors, 2);
238  expect_errors(tano_errors, 2);
239  expect_errors(tano_errors_twice, 2);
240 }
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
T & set(const std::string &name, bool quiet_mode=false)
MeshBase & mesh
static constexpr std::size_t dim
std::unique_ptr< MooseApp > create(const std::string &app_type, const std::string &name, InputParameters parameters, MPI_Comm COMM_WORLD_IN)
void interpolateReconstruct(CellCenteredMapFunctor< T, Map > &output_functor, const Moose::FunctorBase< T > &input_functor, const unsigned int num_int_recs, const bool weight_with_sf, const std::vector< const FaceInfo *> &faces, const Moose::StateArg &time)
Takes an input functor that can be evaluated at faces, typically by linearly interpolating between ad...
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
StateArg currentState()
auto index_range(const T &sizable)
A functor whose evaluation relies on querying a map where the keys are element ids and the values cor...