diff --git a/src/main.cpp b/src/main.cpp
index 67ea80d135ebb922e9108033c92eaaef6f01c5e2..81a4b78b9c9b354e6cff986181e593eec874cfc5 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -141,7 +141,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=0.7;
+ const double tmax=1.5;
double t=0.;
int itermax=std::numeric_limits<int>::max();
@@ -280,7 +280,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);
@@ -304,7 +304,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
<< ": " << rang::fgB::green << error4 << rang::fg::reset << " \n";
- */
+
/*
double error3 = 0.;
error3 = finite_volumes_diffusion.error_L2_E(unknowns);
@@ -329,12 +329,12 @@ int main(int argc, char *argv[])
{ // 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");
+ double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
+ std::ofstream fout("rho800");
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';
+ 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';
}
}
@@ -347,8 +347,8 @@ int main(int argc, char *argv[])
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] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
- fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
+ //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
}
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 470236c9e073d2670ebd4c2af54d74619bbc62e5..1f26bcac1ae9559bfe59af8a5b34be4904490654 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 = -1.;
- double b = 2.;
+ double a = 0.;
+ double b = 1.;
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 407f07b33da97b78cc1eb02fa2cd70ab653102eb..ed4c99543dab6c381ff54b6965b1f5e816f23633 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -197,8 +197,8 @@ private:
m_ur[r]=invAr(r)*br(r);
});
- m_ur[0]=zero;
- m_ur[m_mesh.numberOfNodes()-1]=zero;
+ //m_ur[0]=zero;
+ //m_ur[m_mesh.numberOfNodes()-1]=zero;
// Kidder
@@ -207,11 +207,11 @@ private:
//m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
//R(t) = x*h(t) a la place de x(t)
- /*
+
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;
}
@@ -386,9 +386,9 @@ public:
// Mise a jour de k
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- //kj(j) = xj[j][0];
+ kj(j) = xj[j][0];
//kj(j) = 0.5;
-
+ /*
if (xj[j][0]<0.7) {
kj[j]=0.;
} else {
@@ -398,7 +398,7 @@ public:
kj[j]=0. ;
}
}
-
+ */
});
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 626cc2bf79ac7f6a1f8a73b035156af7dd0f4ddf..8b1206df1c8d6985b77bc260dc1005f8cd3c737e 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -115,7 +115,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)));
@@ -124,7 +124,7 @@ private:
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)));
//m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
-
+ */
// Kidder
@@ -135,11 +135,11 @@ private:
// k = x
//m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][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 ;
}
@@ -178,18 +178,18 @@ 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
//m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
//m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
- /*
+
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 ;
@@ -332,7 +332,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) {
+ 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.;
int l = 0;
@@ -350,8 +351,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)));
+ 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)));
});
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index b1abbde581e6045076a371a33ef3ee6758d27b9e..405e94c1e6933c30dc7a0537f1de24aeb8def18c 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -223,27 +223,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);
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -251,8 +251,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);
@@ -272,9 +272,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;
-
+ /*
if (xj[j][0]<0.7) {
m_kj[j]=0.;
} else {
@@ -284,7 +284,7 @@ public:
m_kj[j]=0. ;
}
}
-
+ */
});
// Conditions aux bords de Dirichlet sur u et k
@@ -292,7 +292,7 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
m_kL[0] = 0.;
- m_kR[0] = 0.;
+ m_kR[0] = 1.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));