Select Git revision
Heat5PointsAlgorithm.hpp
-
PATELA Julie authoredPATELA Julie authored
main.cpp 11.30 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 <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=1.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);
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;
}
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 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;
}