diff --git a/CMakeLists.txt b/CMakeLists.txt
index f363cbb10fa2bdcbce618f0f0e04972eff4a8adc..61d2259bd066be830dabfc0c1d190258549ac153 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -426,6 +426,7 @@ target_link_libraries(
PugsLanguageAST
PugsLanguageModules
PugsMesh
+ PugsUtils
PugsLanguageUtils
kokkos
${PARMETIS_LIBRARIES}
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b1d6068829fb5df8bfa681976a75605ce8a6f469..c04f29ad41df122df252cdc9e7009b49241e8ed2 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -64,7 +64,7 @@ struct GlaceScheme
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<Dimension>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
std::shared_ptr<const MeshType> m_mesh;
@@ -75,8 +75,6 @@ struct GlaceScheme
const FunctionSymbolId& p_id)
: m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
{
- MeshDataType mesh_data(*m_mesh);
-
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::vector<BoundaryConditionHandler> bc_list;
@@ -115,24 +113,25 @@ struct GlaceScheme
}
}
- UnknownsType unknowns(mesh_data);
+ UnknownsType unknowns(*m_mesh);
+
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
- unknowns.rhoj() =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
+ unknowns.rhoj() =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
+ mesh_data.xj());
- unknowns.pj() =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
+ unknowns.pj() =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
- unknowns.uj() =
- InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
- mesh_data
- .xj());
+ unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj());
+ }
unknowns.gammaj().fill(1.4);
AcousticSolver acoustic_solver(m_mesh, bc_list);
- const CellValue<const double>& Vj = mesh_data.Vj();
-
const double tmax = 0.2;
double t = 0;
@@ -152,25 +151,33 @@ struct GlaceScheme
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
block_eos.updateEandCFromRhoP();
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const CellValue<const double> Vj = mesh_data.Vj();
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+
+ parallel_for(
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
+ }
VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
while ((t < tmax) and (iteration < itermax)) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
vtk_writer.write(m_mesh,
{NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", m_mesh->xr()},
+ NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
t);
- double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+ double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
if (t + dt > tmax) {
dt = tmax - t;
}
@@ -190,12 +197,16 @@ struct GlaceScheme
std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
<< rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
<< rang::style::reset << '\n';
- vtk_writer.write(m_mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", m_mesh->xr()},
- NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ vtk_writer.write(m_mesh,
+ {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
+ t, true); // forces last output
+ }
}
};
diff --git a/src/main.cpp b/src/main.cpp
index daf0d4edb97c9939c5b25217649452207a54816f..9b853028133fd9035f67f743c62ee2a0a298202a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,425 +1,21 @@
#include <utils/PugsUtils.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/Mesh.hpp>
-
-#include <scheme/AcousticSolver.hpp>
-#include <scheme/BoundaryCondition.hpp>
-
-#include <output/VTKWriter.hpp>
-
-#include <utils/Exceptions.hpp>
-#include <utils/SignalManager.hpp>
-#include <utils/Timer.hpp>
-
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
-#include <scheme/NamedBoundaryDescriptor.hpp>
-#include <scheme/NumberedBoundaryDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-
-#include <mesh/MeshNodeBoundary.hpp>
-
#include <language/PugsParser.hpp>
-#include <mesh/GmshReader.hpp>
-
+#include <mesh/MeshDataManager.hpp>
#include <mesh/SynchronizerManager.hpp>
-#include <rang.hpp>
-
-#include <iostream>
-
-#include <limits>
-#include <map>
-#include <regex>
-
int
main(int argc, char* argv[])
{
std::string filename = initialize(argc, argv);
- std::regex gmsh_regex("(.*).msh");
- if (not std::regex_match(filename, gmsh_regex)) {
- SynchronizerManager::create();
- MeshDataManager::create();
- parser(filename);
- MeshDataManager::destroy();
- SynchronizerManager::destroy();
- finalize();
- return 0;
- }
-
- std::map<std::string, double> method_cost_map;
-
- MeshDataManager::create();
SynchronizerManager::create();
+ MeshDataManager::create();
- if (filename != "") {
- std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
- Timer gmsh_timer;
- gmsh_timer.reset();
- GmshReader gmsh_reader(filename);
- method_cost_map["Mesh building"] = gmsh_timer.seconds();
-
- std::shared_ptr<const IMesh> p_mesh = gmsh_reader.mesh();
-
- switch (p_mesh->dimension()) {
- case 1: {
- std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
- std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
- for (const auto& sym_boundary_name : sym_boundary_name_list) {
- std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
- std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
- SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
- new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
- bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
- }
-
- using ConnectivityType = Connectivity1D;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<1>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
-
- std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
-
- Timer timer;
- timer.reset();
- MeshDataType mesh_data(*mesh);
-
- std::vector<BoundaryConditionHandler> bc_list;
- {
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case IBoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
- dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_node_list = 0;
- i_ref_node_list < mesh->connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
- const RefNodeList& ref_node_list = mesh->connectivity().refItemList<ItemType::node>(i_ref_node_list);
- const RefId& ref = ref_node_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
- new SymmetryBoundaryCondition<MeshType::Dimension>(
- MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("Unknown BCDescription\n");
- }
- }
- }
- }
-
- UnknownsType unknowns(mesh_data);
-
- unknowns.initializeSod();
-
- AcousticSolver acoustic_solver(mesh, bc_list);
-
- using Rd = TinyVector<MeshType::Dimension>;
-
- const CellValue<const double>& Vj = mesh_data.Vj();
-
- const double tmax = 0.2;
- double t = 0;
-
- int itermax = std::numeric_limits<int>::max();
- int iteration = 0;
-
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- VTKWriter vtk_writer("mesh", 0.01);
-
- while ((t < tmax) and (iteration < itermax)) {
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t);
- double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
- if (t + dt > tmax) {
- dt = tmax - t;
- }
- mesh = acoustic_solver.computeNextStep(dt, unknowns);
-
- block_eos.updatePandCFromRhoE();
-
- t += dt;
- ++iteration;
- }
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green
- << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-
- { // gnuplot output for density
- const CellValue<const Rd>& xj = mesh_data.xj();
- const CellValue<const double>& rhoj = unknowns.rhoj();
- std::ofstream fout("rho");
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << rhoj[j] << '\n';
- }
- }
-
- break;
- }
- case 2: {
- // test case boundary condition description
- std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
- std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
- for (const auto& sym_boundary_name : sym_boundary_name_list) {
- std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
- std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
- SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
- new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
- bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
- }
-
- using ConnectivityType = Connectivity2D;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<2>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
-
- std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
-
- Timer timer;
- timer.reset();
- MeshDataType mesh_data(*mesh);
-
- std::vector<BoundaryConditionHandler> bc_list;
- {
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case IBoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
- dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_face_list = 0;
- i_ref_face_list < mesh->connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
- const RefFaceList& ref_face_list = mesh->connectivity().refItemList<ItemType::face>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
- new SymmetryBoundaryCondition<MeshType::Dimension>(
- MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("Unknown BCDescription\n");
- }
- }
- }
- }
-
- UnknownsType unknowns(mesh_data);
-
- unknowns.initializeSod();
-
- AcousticSolver acoustic_solver(mesh, bc_list);
-
- const CellValue<const double>& Vj = mesh_data.Vj();
-
- const double tmax = 0.2;
- double t = 0;
-
- int itermax = std::numeric_limits<int>::max();
- int iteration = 0;
-
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- VTKWriter vtk_writer("mesh", 0.01);
-
- while ((t < tmax) and (iteration < itermax)) {
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t);
- double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
- if (t + dt > tmax) {
- dt = tmax - t;
- }
- mesh = acoustic_solver.computeNextStep(dt, unknowns);
-
- block_eos.updatePandCFromRhoE();
-
- t += dt;
- ++iteration;
- }
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green
- << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
- break;
- }
- case 3: {
- std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
- std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
- for (const auto& sym_boundary_name : sym_boundary_name_list) {
- std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
- std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
- SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
- new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
- bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
- }
-
- using ConnectivityType = Connectivity3D;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<3>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
-
- std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(gmsh_reader.mesh());
-
- Timer timer;
- timer.reset();
- MeshDataType mesh_data(*mesh);
-
- std::vector<BoundaryConditionHandler> bc_list;
- {
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case IBoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
- dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_face_list = 0;
- i_ref_face_list < mesh->connectivity().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
- const RefFaceList& ref_face_list = mesh->connectivity().refItemList<ItemType::face>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
- new SymmetryBoundaryCondition<MeshType::Dimension>(
- MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("Unknown BCDescription\n");
- }
- }
- }
- }
-
- UnknownsType unknowns(mesh_data);
-
- unknowns.initializeSod();
-
- AcousticSolver acoustic_solver(mesh, bc_list);
-
- const CellValue<const double>& Vj = mesh_data.Vj();
-
- const double tmax = 0.2;
- double t = 0;
-
- int itermax = std::numeric_limits<int>::max();
- int iteration = 0;
-
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- VTKWriter vtk_writer("mesh", 0.01);
-
- while ((t < tmax) and (iteration < itermax)) {
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t);
- double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
- if (t + dt > tmax) {
- dt = tmax - t;
- }
- mesh = acoustic_solver.computeNextStep(dt, unknowns);
- block_eos.updatePandCFromRhoE();
-
- t += dt;
- ++iteration;
- }
- vtk_writer.write(mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh->xr()},
- NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green
- << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
- break;
- }
- }
-
- std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
- << __LINE__ << ")\n";
-
- } else {
- throw NormalError("Connectivity1D defined by number of nodes no more implemented\n");
- }
+ parser(filename);
MeshDataManager::destroy();
SynchronizerManager::destroy();
-
finalize();
-
- std::string::size_type size = 0;
- for (const auto& method_cost : method_cost_map) {
- size = std::max(size, method_cost.first.size());
- }
-
- for (const auto& method_cost : method_cost_map) {
- std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
- << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
- }
-
return 0;
}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 6a525a6aa731baff2b901fe6086a8309cc1ad3e4..86e78f5f7fde25656882535258f6fb3779855d01 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -272,7 +272,8 @@ class MeshData : public IMeshData
void
_checkCellVolume()
{
- auto Vj = this->Vj();
+ Assert(m_Vj);
+ auto& Vj = *m_Vj;
bool is_valid = [&] {
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
@@ -357,23 +358,18 @@ class MeshData : public IMeshData
return *m_Vj;
}
- void
- updateAllData()
- {
- ;
- }
-
+ private:
+ // MeshData **must** be constructed through MeshDataManager
+ friend class MeshDataManager;
MeshData(const MeshType& mesh) : m_mesh(mesh) {}
+ public:
MeshData() = delete;
- MeshData(const MeshData&) = default;
- MeshData(MeshData&&) = default;
+ MeshData(const MeshData&) = delete;
+ MeshData(MeshData&&) = delete;
- ~MeshData()
- {
- ;
- }
+ ~MeshData() {}
};
#endif // MESH_DATA_HPP
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
index ae77ee8255f62447f8d49cd8f9a85a7aec6ddd10..281e7517aa704e6c4352769e88f854b86cbc7366 100644
--- a/src/mesh/MeshDataManager.cpp
+++ b/src/mesh/MeshDataManager.cpp
@@ -25,8 +25,8 @@ MeshDataManager::destroy()
if (m_instance->m_mesh_mesh_data_map.size() > 0) {
std::stringstream error;
error << ": some mesh data is still registered\n";
- for (const auto& i_connectivity_synchronizer : m_instance->m_mesh_mesh_data_map) {
- error << " - mesh data " << rang::fgB::magenta << i_connectivity_synchronizer.first << rang::style::reset << '\n';
+ for (const auto& i_mesh_data : m_instance->m_mesh_mesh_data_map) {
+ error << " - mesh data " << rang::fgB::magenta << i_mesh_data.first << rang::style::reset << '\n';
}
throw UnexpectedError(error.str());
}
@@ -46,11 +46,12 @@ MeshDataManager::getMeshData(const Mesh<Connectivity<Dimension>>& mesh)
{
const IMesh* p_mesh = &mesh;
- if (auto connectivity_synchronizer = m_mesh_mesh_data_map.find(p_mesh);
- connectivity_synchronizer != m_mesh_mesh_data_map.end()) {
- return dynamic_cast<MeshData<Dimension>&>(*connectivity_synchronizer->second);
+ if (auto i_mesh_data = m_mesh_mesh_data_map.find(p_mesh); i_mesh_data != m_mesh_mesh_data_map.end()) {
+ return dynamic_cast<MeshData<Dimension>&>(*i_mesh_data->second);
} else {
- std::shared_ptr mesh_data = std::make_shared<MeshData<Dimension>>(mesh);
+ // **cannot** use make_shared since MeshData constructor is **private**
+ std::shared_ptr<MeshData<Dimension>> mesh_data{new MeshData<Dimension>(mesh)};
+
m_mesh_mesh_data_map[p_mesh] = mesh_data;
return *mesh_data;
}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 971c8fd2458c35d55013a084a4919d7dce026ad6..d0eb4d0b699709e0cddff74035e6b2263b7fffd2 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -32,7 +32,7 @@ class AcousticSolver
constexpr static size_t Dimension = MeshType::Dimension;
using MeshDataType = MeshData<Dimension>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
std::shared_ptr<const MeshType> m_mesh;
const typename MeshType::Connectivity& m_connectivity;
@@ -315,8 +315,6 @@ class AcousticSolver
parallel_for(
m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
- MeshDataManager::instance().getMeshData(*m_mesh).updateAllData();
-
m_mesh = std::make_shared<MeshType>(m_mesh->shared_connectivity(), new_xr);
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*m_mesh).Vj();
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 406bc571c20729c49cfc26b90d436d438ebf3a15..ec7e8851ed98b64971810cfccbf2b587615e8939 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -5,20 +5,15 @@
#include <mesh/ItemValue.hpp>
#include <scheme/BlockPerfectGas.hpp>
-template <typename TMeshData>
+template <typename MeshType>
class FiniteVolumesEulerUnknowns
{
public:
- using MeshDataType = TMeshData;
- using MeshType = typename MeshDataType::MeshType;
-
static constexpr size_t Dimension = MeshType::Dimension;
- using Rd = TinyVector<Dimension>;
- private:
- MeshDataType& m_mesh_data;
- const MeshType& m_mesh;
+ using Rd = TinyVector<Dimension>;
+ private:
CellValue<double> m_rhoj;
CellValue<Rd> m_uj;
CellValue<double> m_Ej;
@@ -139,61 +134,16 @@ class FiniteVolumesEulerUnknowns
return m_inv_mj;
}
- void
- initializeSod()
- {
- const CellValue<const Rd>& xj = m_mesh_data.xj();
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- if (xj[j][0] < 0.5) {
- m_rhoj[j] = 1;
- } else {
- m_rhoj[j] = 0.125;
- }
- });
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- if (xj[j][0] < 0.5) {
- m_pj[j] = 1;
- } else {
- m_pj[j] = 0.1;
- }
- });
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_uj[j] = zero; });
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_gammaj[j] = 1.4; });
-
- BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
- block_eos.updateEandCFromRhoP();
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_Ej[j] = m_ej[j] + 0.5 * (m_uj[j], m_uj[j]); });
-
- const CellValue<const double>& Vj = m_mesh_data.Vj();
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_mj[j] = m_rhoj[j] * Vj[j]; });
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_inv_mj[j] = 1. / m_mj[j]; });
- }
-
- FiniteVolumesEulerUnknowns(MeshDataType& mesh_data)
- : m_mesh_data(mesh_data),
- m_mesh(m_mesh_data.mesh()),
- m_rhoj(m_mesh.connectivity()),
- m_uj(m_mesh.connectivity()),
- m_Ej(m_mesh.connectivity()),
- m_ej(m_mesh.connectivity()),
- m_pj(m_mesh.connectivity()),
- m_gammaj(m_mesh.connectivity()),
- m_cj(m_mesh.connectivity()),
- m_mj(m_mesh.connectivity()),
- m_inv_mj(m_mesh.connectivity())
+ FiniteVolumesEulerUnknowns(const MeshType& mesh)
+ : m_rhoj(mesh.connectivity()),
+ m_uj(mesh.connectivity()),
+ m_Ej(mesh.connectivity()),
+ m_ej(mesh.connectivity()),
+ m_pj(mesh.connectivity()),
+ m_gammaj(mesh.connectivity()),
+ m_cj(mesh.connectivity()),
+ m_mj(mesh.connectivity()),
+ m_inv_mj(mesh.connectivity())
{
;
}