Select Git revision
ModuleRepository.hpp
-
Stéphane Del Pino authored
Plugins are loaded through environment variables. Two environment variables are read: PUGS_PLUGIN and PUGS_PLUGIN_DIR - PUGS_PLUGIN is a string literal that contains the filename of the dynamic library that contains plugins. One can provide multiple filename using semicolumn separators. ex. PUGS_PLUGIN="/path/to/my/libplugin1.so;/anotherpath/to/another/libplugin2.so" - PUGS_PLUGIN_DIR is a string literal that contains directory path where plugins (dynamic libraries) are. One can provide multiple path using semicolumn separators. All the dynamic libraries present at the locations are loaded! ex. PUGS_PLUGIN="/path/to/a/plugin/list/;/anotherpath/to/another/plugin/list/"
Stéphane Del Pino authoredPlugins are loaded through environment variables. Two environment variables are read: PUGS_PLUGIN and PUGS_PLUGIN_DIR - PUGS_PLUGIN is a string literal that contains the filename of the dynamic library that contains plugins. One can provide multiple filename using semicolumn separators. ex. PUGS_PLUGIN="/path/to/my/libplugin1.so;/anotherpath/to/another/libplugin2.so" - PUGS_PLUGIN_DIR is a string literal that contains directory path where plugins (dynamic libraries) are. One can provide multiple path using semicolumn separators. All the dynamic libraries present at the locations are loaded! ex. PUGS_PLUGIN="/path/to/a/plugin/list/;/anotherpath/to/another/plugin/list/"
main.cpp 25.05 KiB
#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;
}