diff --git a/src/main.cpp b/src/main.cpp
index 4713c168ee33421c75cff992ae403c9ee38f607b..c8b67fadecce48eed6e8569cce9aa754dcc7236b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
- const double tmax=1.5;
+ const double tmax=0.8;
double t=0.;
int itermax=std::numeric_limits<int>::max();
@@ -537,7 +537,7 @@ int main(int argc, char *argv[])
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);
@@ -575,8 +575,8 @@ 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";
@@ -636,21 +636,34 @@ int main(int argc, char *argv[])
}
}
- /*
+
{ // 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("S");
+ 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])) << ' ' << S0(j) << '\n';
+ 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';
+ }
+ }
+ }
}
Kokkos::finalize();
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 1f26bcac1ae9559bfe59af8a5b34be4904490654..470236c9e073d2670ebd4c2af54d74619bbc62e5 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -80,8 +80,8 @@ public:
m_x0("x0", 1),
m_xmax("xmax", 1)
{
- double a = 0.;
- double b = 1.;
+ double a = -1.;
+ double b = 2.;
m_x0[0][0] = a;
m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells();
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index a1f4dd88f29a5717377b3f77a428333671b887b2..8232fbecd4a08740cfc347243190cbeda9de8065 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -182,15 +182,15 @@ private:
});
- //m_ur[0]=zero;
- //m_ur[m_mesh.numberOfNodes()-1]=zero;
+ m_ur[0]=zero;
+ m_ur[m_mesh.numberOfNodes()-1]=zero;
// CL Kidder
-
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
-
+ */
return m_ur;
}
@@ -363,12 +363,12 @@ public:
});
// Mise a jour de k
-
+ /*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0];
});
-
+ */
// Mise a jour de nu
/*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index b3720dd661bdce8c7d584fddbddcb0ac80da8749..0147c3e7439eeb3d45fe3f2dd92576371fa7e38c 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -113,7 +113,7 @@ private:
});
// Conditions aux bords
- /*
+
int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(0,0);
m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
@@ -121,8 +121,7 @@ private:
cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
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)));
-
- */
+
// Kidder
@@ -131,9 +130,9 @@ private:
//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];
return m_Fl ;
@@ -175,15 +174,15 @@ private:
// Conditions aux bords
- //m_Gl(0) = Fl(0)*uL(0);
- //m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
+ m_Gl(0) = Fl(0)*uL(0);
+ m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder
-
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0);
m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
-
+ */
return m_Gl ;
@@ -231,13 +230,13 @@ private:
});
// Conditions aux bords
- /*
+
int cell_here = face_cells(0,0);
m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0));
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));
- */
+
// Kiddder
@@ -248,8 +247,8 @@ private:
// nu = 0
- m_Bl(0) = 0.;
- m_Bl(m_mesh.numberOfFaces()-1) = 0.;
+ //m_Bl(0) = 0.;
+ //m_Bl(m_mesh.numberOfFaces()-1) = 0.;
// nu = 0.2
//m_Bl(0) = (0.2*3.*h*x0[0][0])/(50.*h*h*h*h);
@@ -416,7 +415,7 @@ public:
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
// double pi = 4.*std::atan(1.);
- double h = std::sqrt(1. - (t*t)/(50./9.));
+ //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.);
@@ -460,8 +459,8 @@ public:
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
// ajout second membre pour kidder (k = x)
- uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
- Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // nu = 0
+ //uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // nu = 0
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(0.2*3.)/(50.*h*h*h*h)); // nu = 0.2
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(6*xj[j][0]+3.)/(100*h*h*h*h)); // nu = (1+x)*0.5
});
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 595646f1fdccf751e921a590465e74c4caa78565..84ce1bbb347bb69310d052b1b6f343ba09fd6405 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -284,27 +284,27 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- /*
+
if (xj[j][0]<0.5) {
m_rhoj[j]=1.;
} else {
m_rhoj[j]=0.125;
}
- */
+
//Kidder
- m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
+ //m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- /*
+
if (xj[j][0]<0.5) {
m_pj[j]=1;
} else {
m_pj[j]=0.1;
}
- */
+
//Kidder
- m_pj[j] = 2.*std::pow(m_rhoj[j],3);
+ //m_pj[j] = 2.*std::pow(m_rhoj[j],3);
});
double pi = 4.*std::atan(1.);
@@ -313,8 +313,8 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- //m_gammaj[j] = 1.4;
- m_gammaj[j] = 3.;
+ m_gammaj[j] = 1.4;
+ //m_gammaj[j] = 3.;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -334,9 +334,9 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_kj[j] = xj[j][0];
+ //m_kj[j] = xj[j][0];
- //m_kj[j] = 0.5;
+ m_kj[j] = 0.2;
// Sod
@@ -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;
@@ -378,7 +378,7 @@ 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.00014*m_kj[j];
+ m_kj[j] = 0.014*m_kj[j];
*/
});
@@ -406,8 +406,8 @@ public:
// Conditions aux bords de Dirichlet sur T et nu
- m_TL[0] = 2.;
- m_TR[0] = 2.;
+ m_TL[0] = m_ej[0];
+ m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
m_nuL[0] = 0.;
m_nuR[0] = 0.;