From 3b3c19daef00a6d120b7cec4fdef76aa7fbabea1 Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Tue, 7 Jul 2020 19:19:24 +0200 Subject: [PATCH] Begin correction of MeshData use --- CMakeLists.txt | 1 + src/language/modules/SchemeModule.cpp | 67 ++-- src/main.cpp | 410 +--------------------- src/mesh/MeshData.hpp | 9 +- src/scheme/AcousticSolver.hpp | 4 +- src/scheme/FiniteVolumesEulerUnknowns.hpp | 76 +--- 6 files changed, 60 insertions(+), 507 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f363cbb10..61d2259bd 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 b1d606882..86a0f7066 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(*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(*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(*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(*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 daf0d4edb..9b8530281 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 6a525a6aa..b924ae302 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -73,6 +73,7 @@ class MeshData : public IMeshData }); m_xj = std::make_shared<CellValue<const Rd>>(xj); } + std::cout << "- computed xj for mesh " << &m_mesh << '\n'; } PUGS_INLINE @@ -96,6 +97,8 @@ class MeshData : public IMeshData Vj[j] = inv_Dimension * sum_cjr_xr; }); m_Vj = std::make_shared<CellValue<const double>>(Vj); + + std::cout << "- computed Vj for mesh " << &m_mesh << '\n'; } PUGS_INLINE @@ -357,12 +360,6 @@ class MeshData : public IMeshData return *m_Vj; } - void - updateAllData() - { - ; - } - MeshData(const MeshType& mesh) : m_mesh(mesh) {} MeshData() = delete; diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 971c8fd24..d0eb4d0b6 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 406bc571c..ec7e8851e 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()) { ; } -- GitLab