diff --git a/src/main.cpp b/src/main.cpp
index 7440d67736e1cc88935e721f02ec5bb6f897da09..bc1424410fe2f7652da9ba489a860a6502fbf96b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -165,15 +165,16 @@ int main(int argc, char *argv[])
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<std::vector<double>> rho_marche(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();
@@ -251,13 +252,13 @@ int main(int argc, char *argv[])
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);
@@ -270,7 +271,7 @@ int main(int argc, char *argv[])
}
}
diff.close();
- */
+
while((t<tmax) and (iteration<itermax)) {
@@ -322,10 +323,10 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
- /*
+
// ECRITURE DANS UN FICHIER
- if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
+ if ((std::fmod(t,0.001) < 0.0001) or (t == tmax)) {
std::string ligne;
@@ -541,12 +542,12 @@ int main(int argc, char *argv[])
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];
@@ -554,11 +555,12 @@ int main(int argc, char *argv[])
}
tempo[i] = t;
i = i + 1;
-
+ */
}
- std::cout << "i = " << i << std::endl;
+ //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";
@@ -574,7 +576,7 @@ int main(int argc, char *argv[])
fout << ' ' << '\n';
}
*/
-
+ /*
std::ofstream fout("cararho");
fout.precision(15);
for (int j=0; j<mesh.numberOfCells(); ++j) {
@@ -588,7 +590,7 @@ int main(int argc, char *argv[])
fout << ' ' << '\n';
}
}
-
+ */
/*
double error1 = 0.;
@@ -602,7 +604,7 @@ int main(int argc, char *argv[])
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);
@@ -640,7 +642,7 @@ int main(int argc, char *argv[])
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";
@@ -667,6 +669,26 @@ int main(int argc, char *argv[])
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
+
+ { // gnuplot output for pression
+ const Kokkos::View<const Rd*> xj = mesh_data.xj();
+ const Kokkos::View<const double*> pj = unknowns.pj();
+ std::ofstream fout("p");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << pj[j] << '\n';
+ }
+ }
+
+ { // gnuplot output for internal energy
+ const Kokkos::View<const Rd*> xj = mesh_data.xj();
+ const Kokkos::View<const double*> ej = unknowns.ej();
+ std::ofstream fout("e");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << ej[j] << '\n';
+ }
+ }
{ // gnuplot output for vitesse
const Kokkos::View<const Rd*> xj = mesh_data.xj();
@@ -676,12 +698,12 @@ int main(int argc, char *argv[])
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(-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';
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
}
@@ -694,12 +716,12 @@ int main(int argc, char *argv[])
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] << ' ' << (-(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';
+ fout << xj[j][0] << ' ' << Ej[j] << '\n';
}
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 074d0f9572446c866192f578e1d4d724bd2517f9..7dc3959daa79b2566c93ec5bca8dd0c989cd5d0c 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -124,17 +124,18 @@ private:
m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
//m_Fl(m_mesh.numberOfFaces()-1) = -((kR(0)/Vj(m_mesh.numberOfCells()-1) + kj(cell_here)/Vj(cell_here))/(1./Vj(cell_here) + 1./Vj(m_mesh.numberOfCells()-1)))*(1./Vl(m_mesh.numberOfFaces()-1))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+
// Kidder
// k = 0.5
//m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
-
+ /*
// k = x
- //double h = std::sqrt(1. - (t*t)/(50./9.));
- //m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
- //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
-
+ double h = std::sqrt(1. - (t*t)/(50./9.));
+ m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
+ m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
+ */
// k = 0
//m_Fl(0,0) = 0.;
//m_Fl(m_mesh.numberOfFaces()-1,0) = 0.;
@@ -244,6 +245,7 @@ private:
cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
m_Bl(m_mesh.numberOfFaces()-1) = -(nuR(0) + nuj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(Tj(cell_here) - TR(0));
//m_Bl(m_mesh.numberOfFaces()-1) = (nuR(0)/Vj(m_mesh.numberOfCells()-1) + nuj(cell_here)/Vj(cell_here))/(1./Vj(m_mesh.numberOfCells()-1) + 1./Vj(cell_here))*((Tj(cell_here)-TR(0))/Vl(m_mesh.numberOfFaces()-1));
+
// Kiddder
/*
@@ -449,11 +451,11 @@ public:
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
- // double pi = 4.*std::atan(1.);
+ double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (t*t)/(50./9.));
// CL en diffusion pure
- // TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
+ //TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
//TL(0) = 1.+ t;
//TR(0) = 3. + t;
@@ -472,8 +474,8 @@ public:
// Mise a jour de la vitesse et de l'energie totale specifique
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- //const int j = j0+1;
+ Kokkos::parallel_for(m_mesh.numberOfCells()-2, KOKKOS_LAMBDA(const int& j0) {
+ const int j = j0+1;
Rd momentum_fluxes = zero;
double energy_fluxes = 0.;
Rd trich = zero;
@@ -594,9 +596,9 @@ public:
double err_u = 0.;
double exact_u = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
+ //exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
//exact_u = std::sin(pi*xj[j][0])*std::exp(-t); // solution exacte cas test k non constant
- //exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
+ exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
//exact_u = xj[j][0];
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
@@ -616,8 +618,8 @@ public:
double err_E = 0.;
double exact_E = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- exact_E = (-(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.;
- //exact_E = ((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.*t)-1.) + 2.;
+ //exact_E = (-(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.;
+ exact_E = ((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.*t)-1.) + 2.;
//exact_E = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1. + t;
err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
}
@@ -634,7 +636,6 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
- //double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (t*t)/(50./9.));
double exacte = std::sqrt((4.*((xj[0][0]*xj[0][0])/(h*h)) + 100.-(xj[0][0]*xj[0][0])/(h*h))/100.)/h;
double erreur = std::abs(exacte - rhoj[0]);
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 9464275392d87641cc1362a8707f06184f0d60a7..19079bb705b1794c29a83c404cfee47aba10478b 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -336,11 +336,11 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = xj[j][0];
- //m_kj[j] = 0.;
+ m_kj[j] = 0.;
/*
// Sod
-
+
// k non regulier
if (xj[j][0]<0.7) {
@@ -364,7 +364,7 @@ public:
//m_kj[j] = 0.0007;
// Re = 10 000
- //m_kj[j] = 0.00014;
+ m_kj[j] = 0.00014;
// Re = 100 000
//m_kj[j] = 0.000014;
@@ -374,12 +374,12 @@ public:
}
}
*/
-
+ /*
// k regulier
int n = 1;
m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n);
- m_kj[j] = 0.14*m_kj[j];
-
+ m_kj[j] = 0.014*m_kj[j];
+ */
});
@@ -398,17 +398,6 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_nuj(j) = 0.5*(1.+xj[j][0]);
m_nuj(j) = 0.;
- /*
- if (xj[j][0]<0.7) {
- m_nuj[j]=0.;
- } else {
- if (xj[j][0]<0.9){
- m_nuj[j]=0.5;
- } else {
- m_nuj[j]=0. ;
- }
- }
- */
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -470,8 +459,8 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_kj[j] = 2.; // k constant, k = 2
- //m_kj[j] = xj[j][0]; // k non constant, k = x
+ //m_kj[j] = 2.; // k constant, k = 2
+ m_kj[j] = xj[j][0]; // k non constant, k = x
});
// Conditions aux bords de Dirichlet sur u et k
@@ -479,12 +468,12 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
//m_uR[0] = xj[0];
- m_kL[0] = 2.;
- m_kR[0] = 2.;
+ m_kL[0] = 0.;
+ m_kR[0] = 1.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_nuj(j) = 0.;
- //m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
+ //m_nuj(j) = 0.;
+ m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -493,10 +482,12 @@ public:
// Conditions aux bords de Dirichlet sur T et nu
- m_TL[0] = m_ej[0];
- m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
- m_nuL[0] = 0.;
- m_nuR[0] = 0.;
+ //m_TL[0] = m_ej[0];
+ //m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
+ m_TL[0] = 2.;
+ m_TR[0] = 2.;
+ m_nuL[0] = 0.5;
+ m_nuR[0] = 1.;
}
*/