Skip to content
Snippets Groups Projects
Select Git revision
  • 49759b2895b05974c91f534d1e857aee44cdd861
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

main.cpp

Blame
  • main.cpp 10.33 KiB
    #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 <Mesh.hpp>
    #include <BoundaryCondition.hpp>
    #include <AcousticSolver.hpp>
    
    #include <TinyVector.hpp>
    #include <TinyMatrix.hpp>
    
    #include <GmshReader.hpp>
    
    #include <CLI/CLI.hpp>
    #include <cassert>
    #include <limits>
    #include <map>
    
    int main(int argc, char *argv[])
    {
      long unsigned number = 10;
      std::string filename;
      {
        CLI::App app{"Pastis help"};
    
        app.add_option("number,-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();
      // }
    
      if (filename != "") {
        std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
        GmshReader gmsh_reader(filename);
    
        typedef Mesh<Connectivity2D> MeshType;
        typedef MeshData<MeshType> MeshDataType;
        typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
    
        MeshType& mesh = 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) {
    	Connectivity2D::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
    	std::cout << i_boundary
    		  << " -> "
    		  << mesh.connectivity().nodesBoundary(i_boundary).second.extent(0)
    		  << '\n';
    	unsigned int ref = nodes_boundary.first;
    	TinyVector<2> normal(0,0);
    	if ((ref == 3) or (ref == 4)) {
    	  normal = TinyVector<2>(0,1);
    	} else {
    	  normal = TinyVector<2>(1,0);
    	}
    	std::cout << "boundary=" << i_boundary << " ref=" << ref << " n=" << normal << '\n';
    	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));
          }
    
          // 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::ofstream gnuplot("sol.gnu");
        for (size_t j=0; j<mesh.numberOfCells(); ++j) {
          for (int r=0; r<3; ++r) {
    	const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)];
    	gnuplot << x[0] << ' ' << x[1] << '\n';
          }
          const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)];
          gnuplot << x[0] << ' ' << x[1] << "\n\n";
        }
        gnuplot.close();
    
        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();
    
      } else {
        // class for acoustic solver test
        Kokkos::Timer timer;
        timer.reset();
        Connectivity1D connectivity(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;
          }
    
          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';
          }
        }
    
      }
    
      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;
    }