main.cpp

#include <iostream>
#include <Kokkos_Core.hpp>
#include <RevisionInfo.hpp>
#include <rang.hpp>
#include <FPEManager.hpp>
#include <SignalManager.hpp>
#include <ConsoleManager.hpp>

// #include <RawKokkosAcousticSolver.hpp>
// #include <MeshLessAcousticSolver.hpp>
// #include <AcousticSolverClass.hpp>
// #include <AcousticSolverTest.hpp>

#include <Connectivity1D.hpp>
#include <Connectivity2D.hpp>
#include <Connectivity3D.hpp>

#include <Mesh.hpp>
#include <BoundaryCondition.hpp>
#include <AcousticSolver.hpp>

#include <VTKWriter.hpp>

#include <TinyVector.hpp>
#include <TinyMatrix.hpp>

#include <BoundaryConditionDescriptor.hpp>

#include <MeshNodeBoundary.hpp>

#include <GmshReader.hpp>
#include <PastisParser.hpp>

#include <CLI/CLI.hpp>
#include <limits>
#include <map>

int main(int argc, char *argv[])
{
  if (argc == 2) {
    parser(argv[1]);
    return 0;
  }

  long unsigned number = 10;
  std::string filename;
  {
    CLI::App app{"Pastis help"};

    app.add_option("-n,--number", number, "Number of cells");//->required();

    app.add_option("filename,-f,--filename", filename, "gmsh file");//->required();

    int threads=-1;
    app.add_option("--threads", threads, "Number of Kokkos threads")->check(CLI::Range(1,std::numeric_limits<decltype(threads)>::max()));

    std::string colorize="auto";
    app.add_set("--colorize", colorize, {"auto", "yes", "no"}, "Colorize console output", true);

    bool disable_fpe = false;
    app.add_flag("--no-fpe", disable_fpe, "Do not trap floating point exceptions");
    bool disable_signals = false;
    app.add_flag("--no-signal", disable_signals, "Do not catches signals");

    std::string pause_on_error="auto";
    app.add_set("--pause-on-error", pause_on_error, {"auto", "yes", "no"}, "Pause for debugging on unexpected error", true);

    std::atexit([](){std::cout << rang::style::reset;});
    try {
      app.parse(argc, argv);
    } catch (const CLI::ParseError &e) {
      return app.exit(e);
    }

    ConsoleManager::init(colorize);
    FPEManager::init(not disable_fpe);
    SignalManager::setPauseForDebug(pause_on_error);
    SignalManager::init(not disable_signals);
  }

  std::cout << "Code version: "
            << rang::style::bold << RevisionInfo::version() << rang::style::reset << '\n';

  std::cout << "-------------------- "
            << rang::fg::green
            << "git info"
            << rang::fg::reset
            <<" -------------------------"
            << '\n';
  std::cout << "tag:  " << rang::fg::reset
            << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
  std::cout << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
  std::cout << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";

  if (RevisionInfo::gitIsClean()) {
    std::cout << rang::fgB::green << "clean" << rang::fg::reset;
  } else {
    std::cout << rang::fgB::red << "dirty" << rang::fg::reset;
  }
  std::cout << ")\n";
  std::cout << "-------------------------------------------------------\n";

  Kokkos::initialize(argc, argv);
  Kokkos::DefaultExecutionSpace::print_configuration(std::cout);

  std::map<std::string, double> method_cost_map;

  // { // Basic function based acoustic solver
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   RawKokkos::AcousticSolver(number);
  //   method_cost_map["RawKokkos"] = timer.seconds();
  // }

  // { // class for acoustic solver (mesh less)
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   MeshLessAcousticSolver acoustic_solver(number);
  //   method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
  // }

  // { // class for acoustic solver
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   AcousticSolverClass acoustic_solver(number);
  //   method_cost_map["AcousticSolverClass"] = timer.seconds();
  // }

  // { // class for acoustic solver test
  //   Kokkos::Timer timer;
  //   timer.reset();
  //   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";
      Kokkos::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->meshDimension()) {
        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));
          }

          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;
          {
            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: {
                  std::cerr << "Unknown BCDescription\n";
                  std::exit(1);
                }
              }
            }
          }

          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

          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';
            }
          }

          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));
          }

          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;
          {
            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: {
                  std::cerr << "Unknown BCDescription\n";
                  std::exit(1);
                }
              }
            }
          }

          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

          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: {
          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));
          }

          typedef Connectivity3D 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;
          {
            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: {
                  std::cerr << "Unknown BCDescription\n";
                  std::exit(1);
                }
              }
            }
          }

          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

          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;
        }
      }

      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_bc0
            = new PressureBoundaryCondition(1,
                                            std::vector<unsigned int>({static_cast<unsigned int>(0u)}));
        std::shared_ptr<PressureBoundaryCondition> bc0(pres_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;
        }

        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();

  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) {
    std::cout << "* ["
      	      << 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;
}