main.cpp

#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 <cmath>

// #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.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();
    Kokkos::View<double*> kj = unknowns.kj();
    Kokkos::View<double*> nuj = unknowns.nuj();
    Kokkos::View<Rd*> uj = unknowns.uj();
    BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);

    Kokkos::View<double*> nuL = unknowns.nuL();
    Kokkos::View<double*> nuR = unknowns.nuR();
    Kokkos::View<double*> kL = unknowns.kL();
    Kokkos::View<double*> kR = unknowns.kR();
    
    double c = 0.;
    c = finite_volumes_diffusion.conservatif(unknowns);
    
    const Kokkos::View<const Rd*> xj = mesh_data.xj();
    int size = 3000;
    std::vector<std::vector<double>> x(size, std::vector<double>(mesh.numberOfCells()));
    std::vector<double> tempo(size);
    int i = 0;

    /*
    // Ecriture des valeurs initiales dans un fichier
    
    const Kokkos::View<const Rd*> xj = mesh_data.xj();
    const Kokkos::View<const Rd*> xr = mesh.xr();
    
    // rho
    std::ofstream fout1("inter1", std::ios::trunc);
    fout1.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      fout1 << std::fixed << xj[j][0] << ' ' << std::fixed << rhoj[j] << '\n';
    }
    fout1.close();

    // u
    std::ofstream fout2("inter2", std::ios::trunc);
    fout2.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      fout2 << std::fixed << xj[j][0] << ' ' << std::fixed << uj[j][0] << '\n';
    }
    fout2.close();

    // e
    std::ofstream fout3("inter3", std::ios::trunc);
    fout3.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      fout3 << std::fixed << xj[j][0] << ' ' << std::fixed << ej[j] << '\n';
    }
    fout3.close();

    // p
    std::ofstream fout4("inter4", std::ios::trunc);
    fout4.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      fout4 << std::fixed << xj[j][0] << ' ' << std::fixed << pj[j] << '\n';
    }
    fout4.close();

    // S
    std::ofstream fout5("inter5", std::ios::trunc);
    fout5.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      fout5 << std::fixed << xj[j][0] << ' ' << std::fixed << std::log(pj[j]*std::pow(rhoj[j],-gammaj[j])) << '\n';
    }
    fout5.close();

    // derivee de u
    std::ofstream fout6("inter6", std::ios::trunc);
    fout6.precision(15);
    // en 1D : mesh.numberOfNodes() = mesh.numberOfCells() - 1
    for (size_t j=0; j<mesh.numberOfNodes()-2; ++j) {
      fout6 << std::fixed << xj[j][0] << ' ' << std::fixed << (uj[j+1][0]-uj[j][0])/(xr[j+2][0]-xr[j+1][0]) << '\n';
    }
    fout6.close();

    // derivee de p
    std::ofstream fout7("inter7", std::ios::trunc);
    fout7.precision(15);
    for (size_t j=0; j<mesh.numberOfNodes()-2; ++j) {
      fout7 << std::fixed << xj[j][0] << ' ' << std::fixed << (pj[j+1]-pj[j])/(xr[j+2][0]-xr[j+1][0]) << '\n';
    }
    fout7.close();

    // derivee de rho
    std::ofstream fout8("inter8", std::ios::trunc);
    fout8.precision(15);
    for (size_t j=0; j<mesh.numberOfNodes()-2; ++j) {
      fout8 << std::fixed << xj[j][0] << ' ' << std::fixed << (rhoj[j+1]-rhoj[j])/(xr[j+2][0]-xr[j+1][0]) << '\n';
    }
    fout8.close();

    // terme baroclyne
    std::ofstream fout9("inter9", std::ios::trunc);
    fout9.precision(15);
    for (size_t j=0; j<mesh.numberOfNodes()-2; ++j) {
      fout9 << std::fixed << xj[j][0] << ' ' << std::fixed << (((rhoj[j+1]-rhoj[j])/(xr[j+2][0]-xr[j+1][0]))*((pj[j+1]-pj[j])/(xr[j+2][0]-xr[j+1][0])))/(rhoj[j]*rhoj[j]) << '\n';
    }
    fout9.close();

    // Fichier temps 
    std::ofstream tempo("temps");
    tempo.precision(5);
    tempo << std::fixed << t << '\n';
    tempo.close();
    
    // Fichier k indicateur
    std::ofstream diff("diffinter");
    diff.precision(5);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      //diff << std::fixed << xj[j][0] << ' ' << std::fixed << kj[j] << '\n';
      if (kj[j]>0.) {
	diff << std::fixed << xj[j][0] << ' ' << std::fixed << 4. << '\n';
      } else {
	diff << std::fixed << xj[j][0] << ' ' << std::fixed << -0.1 << '\n';
      }
    }
    diff.close();
    */


    while((t<tmax) and (iteration<itermax)) {
     
      /*
      // ETAPE 1 DU SPLITTING - EULER
      
      double dt_euler = 0.9*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.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj);
      double t_diff = t-dt_euler;
      
      if (dt_euler <= dt_diff) {
	dt_diff = dt_euler;
	finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
      } else {
	while (t > t_diff) {
	  dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj);
	  if (t_diff+dt_diff > t) {
	    dt_diff = t-t_diff;
	  }
	  finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
	  t_diff += dt_diff;
	}
      }
      */


      // DIFFUSION PURE
      
      double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
      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 DANS UN FICHIER 
      
      if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
    
      std::string ligne;

      // rho
      std::ifstream fint1("inter1");
      std::ofstream fout1("film_rho", std::ios::trunc);
      fout1.precision(15);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(fint1, ligne);
	fout1 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << rhoj[j] << '\n';
      }
      fint1.close();
      fout1.close();
      
      std::ifstream rint1("film_rho");
      std::ofstream rout1("inter1", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(rint1, ligne);
	rout1 << ligne << '\n';
      }
      rint1.close();
      rout1.close();

      // u
      std::ifstream fint2("inter2");
      std::ofstream fout2("film_u", std::ios::trunc);
      fout2.precision(15);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(fint2, ligne);
	fout2 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << uj[j][0] << '\n';
      }
      fint2.close();
      fout2.close();
      
      std::ifstream rint2("film_u");
      std::ofstream rout2("inter2", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(rint2, ligne);
	rout2 << ligne << '\n';
      }
      rint2.close();
      rout2.close();

      // e
      std::ifstream fint3("inter3");
      std::ofstream fout3("film_e", std::ios::trunc);
      fout3.precision(15);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(fint3, ligne);
	fout3 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << ej[j] << '\n';
      }
      fint3.close();
      fout3.close();
      
      std::ifstream rint3("film_e");
      std::ofstream rout3("inter3", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(rint3, ligne);
	rout3 << ligne << '\n';
      }
      rint3.close();
      rout3.close();

      // p
      std::ifstream fint4("inter4");
      std::ofstream fout4("film_p", std::ios::trunc);
      fout4.precision(15);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(fint4, ligne);
	fout4 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << pj[j] << '\n';
      }
      fint4.close();
      fout4.close();
      
      std::ifstream rint4("film_p");
      std::ofstream rout4("inter4", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(rint4, ligne);
	rout4 << ligne << '\n';
      }
      rint4.close();
      rout4.close();

      // S
      std::ifstream fint5("inter5");
      std::ofstream fout5("film_S", std::ios::trunc);
      fout5.precision(15);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(fint5, ligne);
	fout5 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << std::log(pj[j]*std::pow(rhoj[j],-gammaj[j])) << '\n';
      }
      fint5.close();
      fout5.close();
      
      std::ifstream rint5("film_S");
      std::ofstream rout5("inter5", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(rint5, ligne);
	rout5 << ligne << '\n';
      }
      rint5.close();
      rout5.close();

      // derivee de u
      std::ifstream fint6("inter6");
      std::ofstream fout6("film_du", std::ios::trunc);
      fout6.precision(15);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(fint6, ligne);
	fout6 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << (uj[j+1][0]-uj[j][0])/(xr[j+2][0]-xr[j+1][0]) << '\n';
      }
      fint6.close();
      fout6.close();
      
      std::ifstream rint6("film_du");
      std::ofstream rout6("inter6", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(rint6, ligne);
	rout6 << ligne << '\n';
      }
      rint6.close();
      rout6.close();

      // derivee de p
      std::ifstream fint7("inter7");
      std::ofstream fout7("film_dp", std::ios::trunc);
      fout7.precision(15);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(fint7, ligne);
	fout7 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << (pj[j+1]-pj[j])/(xr[j+2][0]-xr[j+1][0]) << '\n';
      }
      fint7.close();
      fout7.close();
      
      std::ifstream rint7("film_dp");
      std::ofstream rout7("inter7", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(rint7, ligne);
	rout7 << ligne << '\n';
      }
      rint7.close();
      rout7.close();

      // derivee de rho
      std::ifstream fint8("inter8");
      std::ofstream fout8("film_drho", std::ios::trunc);
      fout8.precision(15);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(fint8, ligne);
	fout8 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << (rhoj[j+1]-rhoj[j])/(xr[j+2][0]-xr[j+1][0]) << '\n';
      }
      fint8.close();
      fout8.close();
      
      std::ifstream rint8("film_drho");
      std::ofstream rout8("inter8", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(rint8, ligne);
	rout8 << ligne << '\n';
      }
      rint8.close();
      rout8.close();

      // terme baroclyne
      std::ifstream fint9("inter9");
      std::ofstream fout9("film_baro", std::ios::trunc);
      fout9.precision(15);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(fint9, ligne);
	fout9 << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << (((rhoj[j+1]-rhoj[j])/(xr[j+2][0]-xr[j+1][0]))*((pj[j+1]-pj[j])/(xr[j+2][0]-xr[j+1][0])))/(rhoj[j]*rhoj[j])<< '\n';
      }
      fint9.close();
      fout9.close();
      
      std::ifstream rint9("film_baro");
      std::ofstream rout9("inter9", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfNodes()-2; ++j) {
	getline(rint9, ligne);
	rout9 << ligne << '\n';
      }
      rint9.close();
      rout9.close();
      
      // Fichier temps 
      std::ofstream tempo("temps", std::ios::app);
      tempo.precision(5);
      tempo << std::fixed << t << '\n';
      tempo.close();
      
      // Fichier k indicateur
      std::ifstream diffint("diffinter");
      std::ofstream diffout("k", std::ios::trunc);
      diffout.precision(5);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(diffint, ligne);
	//diffout << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << kj[j] << '\n';
	if (kj[j]>0.) {
	  diffout << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << 4. << '\n';
	} else {
	  diffout << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << -0.1 << '\n';
	}
      }
      diffint.close();
      diffout.close();
      
      std::ifstream riffint("k");
      std::ofstream riffout("diffinter", std::ios::trunc);
      for (size_t j = 0; j<mesh.numberOfCells(); ++j) {
	getline(riffint, ligne);
	riffout << ligne << '\n';
      }
      riffint.close();
      riffout.close();
      
      }
      */
 
      // ENTROPY TEST
      //finite_volumes_diffusion.entropie(unknowns);

      /*
      // STOCKAGE COORDONNES ET TEMPS
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
	x[i][j] = xj[j][0];
      }
      tempo[i] = t;
      i = i + 1;
      */

    }

    std::cout << "i = " << i << std::endl;
    std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";

    // CREATION FICHIER x = f(t)
    /*
    std::ofstream fout("cara");
    fout.precision(15);
    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
      for (size_t k=0; k<i; ++k) {
	fout << tempo[k] << ' ' << x[k][j] << '\n';
      }
      fout << ' ' << '\n';
      fout << ' ' << '\n';
    }
    */
    /*
    std::ofstream fout("cara1");
    fout.precision(15);
    for (int j=0; j<mesh.numberOfCells(); ++j) {
      if (j%10 == 0) {
	for (int k=0; k<i; ++k) {
	  if ( (k%10 == 0) or (k == i-1) ) {
	    fout << tempo[k] << ' ' << x[k][j] << '\n';
	  }
	}
	fout << ' ' << '\n';
	fout << ' ' << '\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_rho(unknowns, tmax);

    std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
	      << ":  " << rang::fgB::green << error2 << rang::fg::reset << " \n";
    */

    double err0 = 0.;
    err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);

    std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
	      << ":  " << rang::fgB::green << err0 << 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 error4 = 0.;
    error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);

    std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
	      << ":  " << rang::fgB::green << error4 << rang::fg::reset << " \n";

    double err1 = 0.;
    err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);

    std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
	      << ":  " << rang::fgB::green << err1 << rang::fg::reset << " \n";

    double error3 = 0.;
    error3 = finite_volumes_diffusion.error_L2_E(unknowns,tmax);

    std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
	      << ":  " << rang::fgB::green << error3 << rang::fg::reset << " \n";
    
    double error5 = 0.;
    error5 = finite_volumes_diffusion.error_Linf_E(unknowns, tmax);

    std::cout << "* " << rang::style::underline << "Erreur L infini E" << rang::style::reset
	      << ":  " << rang::fgB::green << error5 << 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("rho");
     fout.precision(15);
     for (size_t j=0; j<mesh.numberOfCells(); ++j) {
       //fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
       fout << xj[j][0] << ' ' << rhoj[j] << '\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(-tmax) <<'\n'; // cas k non constant
       //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
       //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << xj[j][0] << std::endl;

       //fout << xj[j][0] << ' ' << uj[j][0] << '\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.*tmax)-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] << ' ' << xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + tmax + 1. << std::endl;

       //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();
     const Kokkos::View<const double*> S0 = unknowns.S0();
     std::ofstream fout("S2");
     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';
     }
     }
     */
     /*
     { // gnuplot output for k
     const Kokkos::View<const Rd*> xj   = mesh_data.xj();
     const Kokkos::View<const double*> kj = unknowns.kj();
     std::ofstream fout("k");
     fout.precision(15);
     for (size_t j=0; j<mesh.numberOfCells(); ++j) {
       if (kj[j]>0.) {
	 fout << xj[j][0] << ' ' << 5. << '\n';
       } else {
	 fout << xj[j][0] << ' ' << 0. << '\n';
       }
     }
     }
     */
     /*
     { // gnuplot output for vitesse
     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();
     const Kokkos::View<const Rd*> uj = unknowns.uj();
     std::ofstream fout("vitesse");
     fout.precision(15);
     for (size_t j=0; j<mesh.numberOfCells(); ++j) {
       fout << xj[j][0] << ' ' << uj[j][0] + std::sqrt((gammaj[j]*pj[j])/rhoj[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;
}