Select Git revision
SchemeModule.cpp
-
Stéphane Del Pino authoredStéphane Del Pino authored
SchemeModule.cpp 10.61 KiB
#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");
}
}
}}
));
}