10 #include "gtest/gtest.h" 30 #include "libmesh/elem.h" 31 #include "libmesh/tensor_value.h" 32 #include "libmesh/point.h" 33 #include "libmesh/vector_value.h" 34 #include "libmesh/utility.h" 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());
53 h[i] = 1. / num_elem[i];
57 const auto nx = num_elem[i];
59 auto * factory = &app->getFactory();
60 std::string mesh_type =
"MeshGeneratorMesh";
62 std::shared_ptr<MeshGeneratorMesh>
mesh;
68 app->actionWarehouse().mesh() =
mesh;
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;
78 lm_mesh = mesh_gen->generate();
79 mesh->setMeshBase(std::move(lm_mesh));
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());
91 faces[i] = &all_fi[i];
93 auto & lm_mesh =
mesh->getMesh();
97 for (
auto *
const elem : lm_mesh.active_element_ptr_range())
99 const auto centroid = elem->vertex_average();
101 std::cos(centroid(0)) * std::sin(centroid(1)));
102 u[elem->id()] =
value;
106 up_weller(*
mesh,
"up_weller",
true);
108 up_moukalled(*
mesh,
"up_moukalled",
true);
110 up_tano(*
mesh,
"up_tano",
true);
112 for (
const auto & fi : all_fi)
114 auto moukalled_reconstruct = [&fi](
auto & functor,
auto & container)
119 const Point surface_vector = fi.normal() * fi.faceArea();
120 auto product = (uf * fi.dCN()) * surface_vector;
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();
128 moukalled_reconstruct(u, up_moukalled);
135 Real weller_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())
141 const auto elem_id = elem->id();
145 auto compute_elem_error =
146 [elem_id, current_h, elem, &analytic](
auto & container,
auto & error)
148 auto & current = libmesh_map_find(container, elem_id);
149 const auto diff = analytic - current;
150 error += diff * diff * current_h * current_h;
153 const auto elem_point_eval = container(
159 compute_elem_error(up_weller, weller_error);
160 compute_elem_error(up_moukalled, moukalled_error);
161 compute_elem_error(up_tano, tano_error);
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);
175 for (
auto *
const elem : lm_mesh.active_element_ptr_range())
177 const auto elem_id = elem->id();
181 auto compute_elem_error = [elem_id, current_h, &analytic](
auto & container,
auto & error)
183 auto & current = libmesh_map_find(container, elem_id);
184 const auto diff = analytic - current;
185 error += diff * diff * current_h * current_h;
187 compute_elem_error(up_tano, tano_error);
189 tano_error = std::sqrt(tano_error);
190 tano_errors_twice.push_back(tano_error);
193 *
mesh,
"not_restricted",
true);
199 catch (std::runtime_error & e)
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);
206 *
mesh, {1},
"is_restricted",
true);
212 catch (std::runtime_error & e)
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);
217 std::string(e.what()).find(
218 "that subdomain id is not one of the subdomain ids the functor is restricted to") !=
223 std::for_each(h.begin(), h.end(), [](
Real & h_elem) { h_elem = std::log(h_elem); });
225 auto expect_errors = [&h](
auto & errors_arg,
Real expected_error)
228 errors_arg.begin(), errors_arg.end(), [](
Real & error) { error = std::log(error); });
232 const auto & coeffs = fit.getCoefficients();
233 EXPECT_NEAR(coeffs[1], expected_error, .05);
237 expect_errors(moukalled_errors, 2);
238 expect_errors(tano_errors, 2);
239 expect_errors(tano_errors_twice, 2);
This class constructs a functional expansion using a separate series for each Cartesian dimension...
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
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)
TEST(TestReconstruction, Cartesian)
auto index_range(const T &sizable)
void testReconstruction(const Moose::CoordinateSystemType coord_type)
A functor whose evaluation relies on querying a map where the keys are element ids and the values cor...