diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt index 545278f586dffffe68bc7e293f0ce7752981d5d6..703f19857cfb99c0053aecf1af0984dab371d72e 100644 --- a/src/language/modules/CMakeLists.txt +++ b/src/language/modules/CMakeLists.txt @@ -5,6 +5,7 @@ add_library(PugsLanguageModules MathModule.cpp MeshModule.cpp ModuleRepository.cpp + SchemeModule.cpp VTKModule.cpp ) diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp index 30ad1914a77d9659360e3793e35704109f2f27a8..768661e9999a3dd720455f23daa7c252f555d16f 100644 --- a/src/language/modules/ModuleRepository.cpp +++ b/src/language/modules/ModuleRepository.cpp @@ -3,6 +3,7 @@ #include <language/ast/ASTNode.hpp> #include <language/modules/MathModule.hpp> #include <language/modules/MeshModule.hpp> +#include <language/modules/SchemeModule.hpp> #include <language/modules/VTKModule.hpp> #include <language/utils/BuiltinFunctionEmbedder.hpp> #include <language/utils/SymbolTable.hpp> @@ -20,6 +21,7 @@ ModuleRepository::ModuleRepository() this->_subscribe(std::make_unique<MathModule>()); this->_subscribe(std::make_unique<MeshModule>()); this->_subscribe(std::make_unique<VTKModule>()); + this->_subscribe(std::make_unique<SchemeModule>()); } template <typename NameEmbedderMapT, typename EmbedderTableT> diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bc341bc981b4e10dae14da0d13f5051e4cd7e615 --- /dev/null +++ b/src/language/modules/SchemeModule.cpp @@ -0,0 +1,159 @@ +#include <language/modules/SchemeModule.hpp> + +#include <language/utils/BuiltinFunctionEmbedder.hpp> +#include <mesh/Mesh.hpp> +#include <scheme/AcousticSolver.hpp> + +#include <memory> + +/////////// TEMPORARY + +#include <output/VTKWriter.hpp> +#include <scheme/BoundaryConditionDescriptor.hpp> + +template <size_t Dimension> +struct MeshDimensionAdapter +{ + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + using MeshDataType = MeshData<MeshType>; + using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>; + + const MeshType& m_mesh; + + MeshDimensionAdapter(const IMesh& mesh) : m_mesh{dynamic_cast<const MeshType&>(mesh)} + { + std::vector sym_boundary_name_list = [&]() -> std::vector<std::string> { + if constexpr (Dimension == 3) { + return {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"}; + } else if constexpr (Dimension == 2) { + return {"XMIN", "XMAX", "YMIN", "YMAX"}; + } else { + return {"XMIN", "XMAX"}; + } + }(); + + MeshDataType mesh_data(m_mesh); + std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list; + for (const auto& sym_boundary_name : sym_boundary_name_list) { + std::shared_ptr<BoundaryDescriptor> boudary_descriptor = + std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name)); + SymmetryBoundaryConditionDescriptor* sym_bc_descriptor = + new SymmetryBoundaryConditionDescriptor(boudary_descriptor); + + bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor)); + } + + std::vector<BoundaryConditionHandler> bc_list; + { + constexpr ItemType FaceType = [] { + if constexpr (Dimension > 1) { + return ItemType::face; + } else { + return ItemType::node; + } + }(); + + for (const auto& bc_descriptor : bc_descriptor_list) { + switch (bc_descriptor->type()) { + case BoundaryConditionDescriptor::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 < m_mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { + const auto& ref_face_list = m_mesh.connectivity().template refItemList<FaceType>(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>(m_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<MeshDataType> acoustic_solver(mesh_data, 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_" + std::to_string(Dimension), 0.01); + + while ((t < tmax) and (iteration < itermax)) { + 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); + double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj); + if (t + dt > tmax) { + dt = tmax - t; + } + acoustic_solver.computeNextStep(t, dt, unknowns); + + block_eos.updatePandCFromRhoE(); + + t += dt; + ++iteration; + } + 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 + } +}; + +SchemeModule::SchemeModule() +{ + this->_addBuiltinFunction("glace", std::make_shared<BuiltinFunctionEmbedder<void, std::shared_ptr<IMesh>>>( + std::function<void(std::shared_ptr<IMesh>)>{ + + [](std::shared_ptr<const IMesh> p_mesh) -> void { + switch (p_mesh->dimension()) { + case 1: { + MeshDimensionAdapter<1>{*p_mesh}; + break; + } + case 2: { + MeshDimensionAdapter<2>{*p_mesh}; + break; + } + case 3: { + MeshDimensionAdapter<3>{*p_mesh}; + break; + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } + }} + + )); +} diff --git a/src/language/modules/SchemeModule.hpp b/src/language/modules/SchemeModule.hpp new file mode 100644 index 0000000000000000000000000000000000000000..73d675cda9688134709e92a78aaff8fd82d8abbe --- /dev/null +++ b/src/language/modules/SchemeModule.hpp @@ -0,0 +1,21 @@ +#ifndef SCHEME_MODULE_HPP +#define SCHEME_MODULE_HPP + +#include <language/modules/BuiltinModule.hpp> +#include <utils/PugsMacros.hpp> + +class SchemeModule : public BuiltinModule +{ + public: + std::string_view + name() const final + { + return "scheme"; + } + + SchemeModule(); + + ~SchemeModule() = default; +}; + +#endif // SCHEME_MODULE_HPP