diff --git a/src/main.cpp b/src/main.cpp index ea9a348d3da7b1a23f76873d4b80e8a7875bb707..12fe62ea61829affcfea2b2a1de7aad57fb02041 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -119,134 +119,216 @@ int main(int argc, char *argv[]) // AcousticSolverTest acoustic_solver(number); // method_cost_map["AcousticSolverTest"] = timer.seconds(); // } + try { + if (filename != "") { + std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n"; + GmshReader gmsh_reader(filename); + std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh(); + + switch (p_mesh->meshDimension()) { + case 1: { + typedef Connectivity1D ConnectivityType; + typedef Mesh<ConnectivityType> MeshType; + typedef MeshData<MeshType> MeshDataType; + typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; + + const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh()); + + Kokkos::Timer timer; + timer.reset(); + MeshDataType mesh_data(mesh); + + std::vector<BoundaryConditionHandler> bc_list; + // { // quite dirty! + // for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) { + // ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary); + // unsigned int ref = nodes_boundary.first; + // TinyVector<2> normal(0,0); + // if ((ref == 5) or (ref == 6)) { + // normal = TinyVector<2>(0,1); + // } else { + // normal = TinyVector<2>(1,0); + // } + // const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second; + // std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0)); + // for (size_t r=0; r<nodes_ids.extent(0); ++r) { + // node_boundary_vector[r] = nodes_ids[r]; + // } + // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc + // = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal); + // std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); + // bc_list.push_back(BoundaryConditionHandler(bc)); + // } + // } + + UnknownsType unknowns(mesh_data); + + unknowns.initializeSod(); + + AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list); + + typedef TinyVector<MeshType::dimension> Rd; + + const Kokkos::View<const double*> Vj = mesh_data.Vj(); + const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); + + const double tmax=0.2; + double t=0; + + int itermax=std::numeric_limits<int>::max(); + int iteration=0; + + Kokkos::View<double*> rhoj = unknowns.rhoj(); + Kokkos::View<double*> ej = unknowns.ej(); + Kokkos::View<double*> pj = unknowns.pj(); + Kokkos::View<double*> gammaj = unknowns.gammaj(); + Kokkos::View<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, 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, t, true); // forces last output - if (filename != "") { - std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n"; - GmshReader gmsh_reader(filename); - std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh(); + std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset + << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; - switch (p_mesh->meshDimension()) { - case 1: { - typedef Connectivity1D ConnectivityType; - typedef Mesh<ConnectivityType> MeshType; - typedef MeshData<MeshType> MeshDataType; - typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; + method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); - const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh()); + break; + } + case 2: { + typedef Connectivity2D ConnectivityType; + typedef Mesh<ConnectivityType> MeshType; + typedef MeshData<MeshType> MeshDataType; + typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; + + const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh()); + + Kokkos::Timer timer; + timer.reset(); + MeshDataType mesh_data(mesh); + + std::vector<BoundaryConditionHandler> bc_list; + { // quite dirty! + for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) { + ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary); + const RefId& ref = nodes_boundary.first; + TinyVector<2> normal(0,0); + if ((ref.tagName()== std::string("XMIN")) or (ref.tagName()=="XMAX")) { + normal = TinyVector<2>(1,0); + } else if ((ref.tagName()=="YMIN") or (ref.tagName()=="YMAX")) { + normal = TinyVector<2>(0,1); + } else { + std::cout << rang::fg::red << "ignoring " << ref << rang::style::reset << '\n'; + break; + } + + const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second; + std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0)); + for (size_t r=0; r<nodes_ids.extent(0); ++r) { + node_boundary_vector[r] = nodes_ids[r]; + } + SymmetryBoundaryCondition<MeshType::dimension>* sym_bc + = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal); + std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); + bc_list.push_back(BoundaryConditionHandler(bc)); + } + } - Kokkos::Timer timer; - timer.reset(); - MeshDataType mesh_data(mesh); + UnknownsType unknowns(mesh_data); - std::vector<BoundaryConditionHandler> bc_list; - // { // quite dirty! - // for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) { - // ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary); - // unsigned int ref = nodes_boundary.first; - // TinyVector<2> normal(0,0); - // if ((ref == 5) or (ref == 6)) { - // normal = TinyVector<2>(0,1); - // } else { - // normal = TinyVector<2>(1,0); - // } - // const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second; - // std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0)); - // for (size_t r=0; r<nodes_ids.extent(0); ++r) { - // node_boundary_vector[r] = nodes_ids[r]; - // } - // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc - // = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal); - // std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); - // bc_list.push_back(BoundaryConditionHandler(bc)); - // } - // } + unknowns.initializeSod(); - UnknownsType unknowns(mesh_data); + AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list); - unknowns.initializeSod(); + typedef TinyVector<MeshType::dimension> Rd; - AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list); + const Kokkos::View<const double*> Vj = mesh_data.Vj(); + const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); - typedef TinyVector<MeshType::dimension> Rd; + const double tmax=0.2; + double t=0; - const Kokkos::View<const double*> Vj = mesh_data.Vj(); - const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); + int itermax=std::numeric_limits<int>::max(); + int iteration=0; - const double tmax=0.2; - double t=0; + Kokkos::View<double*> rhoj = unknowns.rhoj(); + Kokkos::View<double*> ej = unknowns.ej(); + Kokkos::View<double*> pj = unknowns.pj(); + Kokkos::View<double*> gammaj = unknowns.gammaj(); + Kokkos::View<double*> cj = unknowns.cj(); - int itermax=std::numeric_limits<int>::max(); - int iteration=0; + BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - Kokkos::View<double*> rhoj = unknowns.rhoj(); - Kokkos::View<double*> ej = unknowns.ej(); - Kokkos::View<double*> pj = unknowns.pj(); - Kokkos::View<double*> gammaj = unknowns.gammaj(); - Kokkos::View<double*> cj = unknowns.cj(); + VTKWriter vtk_writer("mesh", 0.01); - BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); + while((t<tmax) and (iteration<itermax)) { + vtk_writer.write(mesh, t); + double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj); + if (t+dt>tmax) { + dt=tmax-t; + } + acoustic_solver.computeNextStep(t,dt, unknowns); - VTKWriter vtk_writer("mesh", 0.01); + block_eos.updatePandCFromRhoE(); - while((t<tmax) and (iteration<itermax)) { - vtk_writer.write(mesh, t); - double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj); - if (t+dt>tmax) { - dt=tmax-t; - } - acoustic_solver.computeNextStep(t,dt, unknowns); + t += dt; + ++iteration; + } + vtk_writer.write(mesh, t, true); // forces last output - block_eos.updatePandCFromRhoE(); + std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset + << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; - t += dt; - ++iteration; + method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); + break; + } + case 3: { + break; + } } - vtk_writer.write(mesh, 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"; + std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n"; - method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); - - break; - } - case 2: { - typedef Connectivity2D ConnectivityType; - typedef Mesh<ConnectivityType> MeshType; + } else { + // class for acoustic solver test + Kokkos::Timer timer; + timer.reset(); + std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number)); + typedef Mesh<Connectivity1D> MeshType; typedef MeshData<MeshType> MeshDataType; typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; - const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh()); - - Kokkos::Timer timer; - timer.reset(); + MeshType mesh(connectivity); MeshDataType mesh_data(mesh); std::vector<BoundaryConditionHandler> bc_list; { // quite dirty! - for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) { - ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary); - const RefId& ref = nodes_boundary.first; - TinyVector<2> normal(0,0); - if ((ref.tagName()== std::string("XMIN")) or (ref.tagName()=="XMAX")) { - normal = TinyVector<2>(1,0); - } else if ((ref.tagName()=="YMIN") or (ref.tagName()=="YMAX")) { - normal = TinyVector<2>(0,1); - } else { - std::cout << rang::fg::red << "ignoring " << ref << rang::style::reset << '\n'; - break; - } - - const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second; - std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0)); - for (size_t r=0; r<nodes_ids.extent(0); ++r) { - node_boundary_vector[r] = nodes_ids[r]; - } - SymmetryBoundaryCondition<MeshType::dimension>* sym_bc - = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal); - std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); - bc_list.push_back(BoundaryConditionHandler(bc)); - } + SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0 + = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}), + TinyVector<1>(-1)); + std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0); + bc_list.push_back(BoundaryConditionHandler(bc0)); + + PressureBoundaryCondition* pres_bc1 + = new PressureBoundaryCondition(0.1, + std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())})); + std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1); + bc_list.push_back(BoundaryConditionHandler(bc1)); } UnknownsType unknowns(mesh_data); @@ -274,14 +356,12 @@ int main(int argc, char *argv[]) BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - VTKWriter vtk_writer("mesh", 0.01); - while((t<tmax) and (iteration<itermax)) { - vtk_writer.write(mesh, 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(); @@ -289,101 +369,26 @@ int main(int argc, char *argv[]) t += dt; ++iteration; } - vtk_writer.write(mesh, 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: { - break; - } - } - - std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n"; - - } else { - // class for acoustic solver test - Kokkos::Timer timer; - timer.reset(); - std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number)); - typedef Mesh<Connectivity1D> MeshType; - typedef MeshData<MeshType> MeshDataType; - typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; - - MeshType mesh(connectivity); - MeshDataType mesh_data(mesh); - - std::vector<BoundaryConditionHandler> bc_list; - { // quite dirty! - SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0 - = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}), - TinyVector<1>(-1)); - std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0); - bc_list.push_back(BoundaryConditionHandler(bc0)); - - PressureBoundaryCondition* pres_bc1 - = new PressureBoundaryCondition(0.1, - std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())})); - std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1); - bc_list.push_back(BoundaryConditionHandler(bc1)); - } - UnknownsType unknowns(mesh_data); - - unknowns.initializeSod(); - - AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list); - - typedef TinyVector<MeshType::dimension> Rd; - - const Kokkos::View<const double*> Vj = mesh_data.Vj(); - const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); - - const double tmax=0.2; - double t=0; - - int itermax=std::numeric_limits<int>::max(); - int iteration=0; - - Kokkos::View<double*> rhoj = unknowns.rhoj(); - Kokkos::View<double*> ej = unknowns.ej(); - Kokkos::View<double*> pj = unknowns.pj(); - Kokkos::View<double*> gammaj = unknowns.gammaj(); - Kokkos::View<double*> cj = unknowns.cj(); - - BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - - while((t<tmax) and (iteration<itermax)) { - double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj); - if (t+dt>tmax) { - dt=tmax-t; + { // gnuplot output for density + const Kokkos::View<const Rd*> xj = mesh_data.xj(); + const Kokkos::View<const double*> rhoj = unknowns.rhoj(); + std::ofstream fout("rho"); + for (size_t j=0; j<mesh.numberOfCells(); ++j) { + fout << xj[j][0] << ' ' << rhoj[j] << '\n'; + } } - acoustic_solver.computeNextStep(t,dt, unknowns); - - block_eos.updatePandCFromRhoE(); - - t += dt; - ++iteration; } - - 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 Kokkos::View<const Rd*> xj = mesh_data.xj(); - const Kokkos::View<const double*> rhoj = unknowns.rhoj(); - std::ofstream fout("rho"); - for (size_t j=0; j<mesh.numberOfCells(); ++j) { - fout << xj[j][0] << ' ' << rhoj[j] << '\n'; - } - } - + } + catch (const AssertError& error) { + std::cerr << error << '\n'; + std::exit(1); } Kokkos::finalize(); diff --git a/src/utils/PastisAssert.hpp b/src/utils/PastisAssert.hpp index a6342a41ee604a861f928d6ede2e9592df0f6bae..ea92e3175fbbb1c28763b684ea09167c18489cb3 100644 --- a/src/utils/PastisAssert.hpp +++ b/src/utils/PastisAssert.hpp @@ -1,14 +1,6 @@ #ifndef PASTIS_ASSERT_HPP #define PASTIS_ASSERT_HPP -#ifdef NDEBUG - -#warning check that this test does not degrade performaces -#define Assert(assertion) \ - if (not (assertion)) false; - -#else // NDEBUG - #include <rang.hpp> #include <iostream> #include <string> @@ -53,6 +45,14 @@ class AssertError ~AssertError() = default; }; +#ifdef NDEBUG + +#warning check that this test does not degrade performaces +#define Assert(assertion) \ + if (not (assertion)) false; + +#else // NDEBUG + #define Assert(assertion) \ if (not (assertion)) { \ throw AssertError(__FILE__, \