#include <language/modules/SchemeModule.hpp> #include <language/utils/BuiltinFunctionEmbedder.hpp> #include <language/utils/TypeDescriptor.hpp> #include <mesh/Mesh.hpp> #include <scheme/AcousticSolver.hpp> #include <scheme/IBoundaryConditionDescriptor.hpp> #include <scheme/IBoundaryDescriptor.hpp> #include <scheme/NamedBoundaryDescriptor.hpp> #include <scheme/NumberedBoundaryDescriptor.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp> #include <memory> /////////// TEMPORARY #include <output/VTKWriter.hpp> template <size_t Dimension> struct GlaceScheme { using ConnectivityType = Connectivity<Dimension>; using MeshType = Mesh<ConnectivityType>; using MeshDataType = MeshData<MeshType>; using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>; const MeshType& m_mesh; GlaceScheme(const IMesh& mesh, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_0, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_1, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_2, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_3, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_4, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_5) : m_mesh{dynamic_cast<const MeshType&>(mesh)} { MeshDataType mesh_data(m_mesh); std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>> bc_descriptor_list; if constexpr (Dimension == 1) { bc_descriptor_list = {boundary_0, boundary_1}; } else if constexpr (Dimension == 2) { bc_descriptor_list = {boundary_0, boundary_1, boundary_2, boundary_3}; } else if constexpr (Dimension == 3) { bc_descriptor_list = {boundary_0, boundary_1, boundary_2, boundary_3, boundary_4, boundary_5}; } std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n'; 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 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 < 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->_addTypeDescriptor( std::make_shared<TypeDescriptor>(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>.typeName())); this->_addTypeDescriptor(std::make_shared<TypeDescriptor>( ast_node_data_type_from<std::shared_ptr<const IBoundaryConditionDescriptor>>.typeName())); this->_addBuiltinFunction("boundaryName", std::make_shared< BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>, std::string>>( std::function<std::shared_ptr<const IBoundaryDescriptor>(std::string)>{ [](const std::string& boundary_name) -> std::shared_ptr<const IBoundaryDescriptor> { return std::make_shared<NamedBoundaryDescriptor>(boundary_name); }} )); this->_addBuiltinFunction("boundaryTag", std::make_shared< BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>, int64_t>>( std::function<std::shared_ptr<const IBoundaryDescriptor>(int64_t)>{ [](int64_t boundary_tag) -> std::shared_ptr<const IBoundaryDescriptor> { return std::make_shared<NumberedBoundaryDescriptor>(boundary_tag); }} )); this ->_addBuiltinFunction("symmetry", std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryDescriptor>>>( std::function<std::shared_ptr<const IBoundaryConditionDescriptor>( std::shared_ptr<const IBoundaryDescriptor>)>{ [](std::shared_ptr<const IBoundaryDescriptor> boundary) -> std::shared_ptr<const IBoundaryConditionDescriptor> { return std::make_shared<SymmetryBoundaryConditionDescriptor>(boundary); }} )); this->_addBuiltinFunction("glace", std::make_shared<BuiltinFunctionEmbedder< void, std::shared_ptr<const IMesh>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>>>( std::function<void(std::shared_ptr<const IMesh>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>, std::shared_ptr<const IBoundaryConditionDescriptor>)>{ [](std::shared_ptr<const IMesh> p_mesh, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_0, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_1, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_2, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_3, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_4, std::shared_ptr<const IBoundaryConditionDescriptor> boundary_5) -> void { switch (p_mesh->dimension()) { case 1: { GlaceScheme<1>{*p_mesh, boundary_0, boundary_1, boundary_2, boundary_3, boundary_4, boundary_5}; break; } case 2: { GlaceScheme<2>{*p_mesh, boundary_0, boundary_1, boundary_2, boundary_3, boundary_4, boundary_5}; break; } case 3: { GlaceScheme<3>{*p_mesh, boundary_0, boundary_1, boundary_2, boundary_3, boundary_4, boundary_5}; break; } default: { throw UnexpectedError("invalid mesh dimension"); } } }} )); }