Select Git revision
ASTSymbolTableBuilder.cpp
-
Stéphane Del Pino authored
This includes an important change in the code since one cannot reuse a variable name in a nested block (this means that if a variable is declared, one cannot declare a variable with a same name even in a nested block, while the first variable exists).
Stéphane Del Pino authoredThis includes an important change in the code since one cannot reuse a variable name in a nested block (this means that if a variable is declared, one cannot declare a variable with a same name even in a nested block, while the first variable exists).
main.cpp 16.03 KiB
#include <PastisUtils.hpp>
#include <PastisOStream.hpp>
#include <rang.hpp>
#include <Connectivity.hpp>
#include <Mesh.hpp>
#include <BoundaryCondition.hpp>
#include <AcousticSolver.hpp>
#include <VTKWriter.hpp>
#include <Timer.hpp>
#include <TinyVector.hpp>
#include <TinyMatrix.hpp>
#include <BoundaryConditionDescriptor.hpp>
#include <MeshNodeBoundary.hpp>
#include <GmshReader.hpp>
#include <PastisParser.hpp>
#include <SynchronizerManager.hpp>
#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)) {
parser(filename);
return 0;
}
std::map<std::string, double> method_cost_map;
SynchronizerManager::create();
if (filename != "") {
pout() << "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<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<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));
}
using ConnectivityType = Connectivity1D;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<MeshType>; using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_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 BoundaryConditionDescriptor::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().numberOfRefNodeList();
++i_ref_node_list) {
const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(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: {
perr() << "Unknown BCDescription\n";
std::exit(1);
}
}
}
}
UnknownsType unknowns(mesh_data);
unknowns.initializeSod();
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, 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;
}
acoustic_solver.computeNextStep(t,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
pout() << "* " << 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<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));
}
using ConnectivityType = Connectivity2D;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<MeshType>;
using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_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 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<mesh.connectivity().numberOfRefFaceList();
++i_ref_face_list) {
const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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: {
perr() << "Unknown BCDescription\n";
std::exit(1);
}
}
}
}
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", 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;
}
acoustic_solver.computeNextStep(t,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
pout() << "* " << 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<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));
}
using ConnectivityType = Connectivity3D;
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<MeshType>;
using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_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 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<mesh.connectivity().numberOfRefFaceList();
++i_ref_face_list) {
const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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: {
perr() << "Unknown BCDescription\n";
std::exit(1);
}
}
}
}
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", 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;
}
acoustic_solver.computeNextStep(t,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
pout() << "* " << 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;
}
}
pout() << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
} else {
perr() << "Connectivity1D defined by number of nodes no more implemented\n";
std::exit(0);
}
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) {
pout() << "* ["
<< 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;
}