#include <iostream> #include <fstream> #include <iomanip> #include <Kokkos_Core.hpp> #include <RevisionInfo.hpp> #include <rang.hpp> #include <FPEManager.hpp> #include <SignalManager.hpp> #include <ConsoleManager.hpp> #include <string> // #include <RawKokkosAcousticSolver.hpp> // #include <MeshLessAcousticSolver.hpp> // #include <AcousticSolverClass.hpp> // #include <AcousticSolverTest.hpp> #include <Connectivity1D.hpp> #include <Mesh.hpp> #include <AcousticSolver.hpp> #include <FiniteVolumesDiffusion.hpp> #include <TinyVector.hpp> #include <TinyMatrix.hpp> #include <CLI/CLI.hpp> #include <cassert> #include <limits> #include <map> int main(int argc, char *argv[]) { long unsigned number = 10; { CLI::App app{"Pastis help"}; app.add_option("number,-n,--number", number, "Number of cells");//->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(); // } { 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); UnknownsType unknowns(mesh_data); unknowns.initializeSod(); AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns); FiniteVolumesDiffusion<MeshDataType> finite_volumes_diffusion(mesh_data, unknowns); 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.5; 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(); Kokkos::View<double*> kj = unknowns.kj(); BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); double c = 0.; c = finite_volumes_diffusion.conservatif(unknowns); // Ecriture des valeurs initiales de rho dans un fichier const Kokkos::View<const Rd*> xj = mesh_data.xj(); std::ofstream fout("film_rho", std::ios::trunc); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << std::fixed << xj[j][0] << ' ' << std::fixed << rhoj[j] << '\n'; } while((t<tmax) and (iteration<itermax)) { // ETAPE 1 DU SPLITTING - EULER double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); if (t+dt_euler > tmax) { dt_euler = tmax-t; } acoustic_solver.computeNextStep(t,dt_euler, unknowns); t += dt_euler; // ETAPE 2 DU SPLITTING - DIFFUSION double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (dt_euler <= dt_diff) { dt_diff = dt_euler; finite_volumes_diffusion.computeNextStep(t, dt_diff, unknowns); } else { double t_diff = t-dt_euler; while (t > t_diff) { finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (t_diff+dt_diff > t) { dt_diff = t-t_diff; } t_diff += dt_diff; } } // DIFFUSION PURE /* double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj); if (t+dt > tmax) { dt = tmax-t; } finite_volumes_diffusion.computeNextStep(t, dt, unknowns); t += dt; */ block_eos.updatePandCFromRhoE(); ++iteration; std::cout << "temps t : " << t << std::endl; // Ecriture des valeurs de rho par pas de temps dans le fichier std::ofstream fout("film_rho", std::ios::app); fout.precision(15); fout.clear(); fout.seekp(0, std::ios::beg); //int const position = 0; for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout.seekp(35*iteration*j); //position = fout.tellp(); //std::cout << "Nous sommes au " << fout.tellp() << "eme caractere du fichier" << std::endl; fout << ' ' << std::fixed << rhoj[j] << '\n'; //std::cin.ignore(35*iteration, '\n'); } } std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; /* double error1 = 0.; error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); std::cout << "* " << rang::style::underline << "Erreur L2 rho" << rang::style::reset << ": " << rang::fgB::green << error1 << rang::fg::reset << " \n"; double error2 = 0.; error2 = finite_volumes_diffusion.error_Linf(unknowns, tmax); std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset << ": " << rang::fgB::green << error2 << rang::fg::reset << " \n"; double error = 0.; error = finite_volumes_diffusion.error_L2_u(unknowns, tmax); std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset << ": " << rang::fgB::green << error << rang::fg::reset << " \n"; */ /* double error3 = 0.; error3 = finite_volumes_diffusion.error_L2_E(unknowns); std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset << ": " << rang::fgB::green << error3 << rang::fg::reset << " \n"; */ std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset << ": " << rang::fgB::green << c << rang::fg::reset << " \n"; double cons = 0.; cons = finite_volumes_diffusion.conservatif(unknowns); std::cout << "* " << rang::style::underline << "Resultat conservativite rho E" << rang::style::reset << ": " << rang::fgB::green << cons << rang::fg::reset << " \n"; //method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds(); { // gnuplot output for density const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> rhoj = unknowns.rhoj(); //double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); std::ofstream fout("rho1.5"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << xj[j][0] << ' ' << rhoj[j] << '\n'; // Kidder // fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; } } { // gnuplot output for vitesse const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> uj = unknowns.uj(); //double pi = 4.*std::atan(1.); std::ofstream fout("u"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2) <<'\n'; //cas k constant //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant fout << xj[j][0] << ' ' << uj[j][0] << '\n'; //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; } } { // gnuplot output for energy const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> Ej = unknowns.Ej(); //double pi = 4.*std::atan(1.); //double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); std::ofstream fout("E"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (-(std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))+(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-1.) + 2. <<'\n'; // cas k constant //fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2. <<'\n' ; // cas k non constant //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder fout << xj[j][0] << ' ' << Ej[j] << '\n'; } } { // gnuplot output for entropy (gaz parfait en prenant cv = 1)) const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> rhoj = unknowns.rhoj(); const Kokkos::View<const double*> pj = unknowns.pj(); const Kokkos::View<const double*> gammaj = unknowns.gammaj(); std::ofstream fout("S"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << xj[j][0] << ' ' << std::log(pj[j]*std::pow(rhoj[j],-gammaj[j])) << '\n'; } } } Kokkos::finalize(); std::cout << "----------------------\n"; 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; }