diff --git a/src/main.cpp b/src/main.cpp
index f31a11b9a5c4319a70a571ba6072f231dde0bf8f..b6afd5ddbe18c22a0f65d2ab7d5ff822db454d1d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -338,7 +338,11 @@ int main(int argc, char *argv[])
// NAVIER-STOKES SANS SPLITTING
- double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+ double dt = 0.1*acoustic_solver.acoustic_dt(Vj, cj); // pour le cas xi = 0
+ double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+ if (dt_diff < dt) {
+ dt = dt_diff;
+ }
if (t+dt > tmax) {
dt = tmax-t;
}
@@ -356,6 +360,7 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
+ // std::exit(0);
// Animations
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index d5a608a3a2133e47a2cfed91e743d9d228b3b06c..383c562b09cc2a6432f5de0e6f26bc852ee84e3f 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -82,8 +82,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/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 60f4ddb6baa20d842c42df4665417cb7d3b9ef0b..a3b67c8c4a29f317b54bf8f42d59c91bdc092348 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -299,28 +299,28 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// TAC
- /*
+
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){
// TAC
- /*
+
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.);
@@ -330,9 +330,9 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// TAC
- //m_gammaj[j] = 1.4;
+ m_gammaj[j] = 1.4;
// Kidder
- m_gammaj[j] = 3.;
+ //m_gammaj[j] = 3.;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -353,7 +353,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// Differents k (xi)
//m_kj[j] = xj[j][0];
- m_kj[j] = 0.5;
+ m_kj[j] = 0.014;
// TAC
@@ -408,8 +408,8 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
- m_kL[0] = 0.5;
- m_kR[0] = 0.5;
+ m_kL[0] = 0.014;
+ m_kR[0] = 0.014;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index 646eb3ea5e7ad2c9784f6ea134d701a8ddb8369b..cacf6f5f6e76630716f16654899ba4ed3ae14b8b 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -78,6 +78,7 @@ private:
{
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
m_rhocj[j] = rhoj[j]*cj[j];
+ //m_rhocj[j] = 1.;
});
return m_rhocj;
}
@@ -159,20 +160,22 @@ private:
m_ur[r]=invAr(r)*br(r);
});
+ // --- CL ---
- //m_ur[0]=zero;
- //m_ur[m_mesh.numberOfNodes()-1]=zero;
+ m_ur[0]=zero;
+ m_ur[m_mesh.numberOfNodes()-1]=zero;
//m_ur[0] = x0;
//m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// 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;
}
@@ -315,7 +318,10 @@ public:
} else {
PTj(j) = pj(j) - kj(j)*2.*sum/sum1;
}
+
+ });
*/
+
/*
// Calcul de PT (2eme essai, symetrisation)
const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
@@ -324,22 +330,7 @@ public:
const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces();
const Kokkos::View<const unsigned short*> cell_nb_faces
= m_connectivity.cellNbFaces();
-
- // gnuplot output for vitesse
- std::ofstream fout("uj");
- fout.precision(15);
- for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << uj[j][0] << '\n';
- }
- // gnuplot output for vitesse riemann
- std::ofstream fout1("ur");
- fout1.precision(15);
- for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
- fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
- }
-
-
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
std::vector<double> stock(2);
@@ -355,43 +346,43 @@ public:
stock[l] = sum/sum2;
}
if (j == 0) {
- //PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
- PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+ PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
+ //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
} else if (j == m_mesh.numberOfCells()-1) {
- PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
- //PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
+ //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+ PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
} else {
PTj(j) = pj(j) - kj(j)*(stock[1]-stock[0])/Vj(j);
}
+
+ });
*/
- // Calcul de PT (3eme essai, avec uR du solveur de Riemann)
- computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, t);
-
- // gnuplot output for vitesse
- std::ofstream fout("uj");
- fout.precision(15);
+ std::ofstream fout2("pj");
+ fout2.precision(15);
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ fout2 << xj[j][0] << ' ' << pj[j] << '\n';
}
-
- // gnuplot output for vitesse riemann
- std::ofstream fout1("ur");
- fout1.precision(15);
- for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
- fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
+ std::ofstream fout3("pTj");
+ fout3.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+ fout3 << xj[j][0] << ' ' << PTj[j] << '\n';
}
+ // Calcul de PT (3eme essai, avec uR du solveur de Riemann)
+ computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, t);
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
if (j == 0) {
- //PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
- PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+ PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
+ //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+ //PTj(j) = pj(j) - kj(j)*(ur[cell_nodes(j,1)][0])/Vj(j);
} else if (j == m_mesh.numberOfCells()-1) {
- PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
- //PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
+ PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
+ //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+ //PTj(j) = pj(j) - kj(j)*(-ur[cell_nodes(j,0)][0])/Vj(j);
} else {
-
double sum = 0;
for (int k=0; k<cell_nb_nodes(j); ++k) {
int node_here = cell_nodes(j,k);
@@ -400,10 +391,10 @@ public:
PTj(j) = pj(j) - kj(j)*sum/Vj(j);
//std::cout << PTj(j) << std::endl;
//PTj(j) = pj(j) - kj(j)*(ur[cell_nodes(j,0)][0]-ur[cell_nodes(j,1)][0])/Vj(j);
- }
-
+ }
+
});
-
+
// Calcule les flux
computeExplicitFluxes(xr, xj, rhoj, uj, PTj, cj, Vj, Cjr, t);
@@ -422,7 +413,7 @@ public:
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
// ajout second membre pour kidder (k cst)
- Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*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)));
@@ -447,13 +438,28 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
rhoj[j] = mj[j]/Vj[j];
});
-
+
+ // gnuplot output for vitesse
+ std::ofstream fout("uj");
+ fout.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ }
+
+ // gnuplot output for vitesse riemann
+ std::ofstream fout1("ur");
+ fout1.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
+ }
+
// Mise a jour de k
/*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0];
});
*/
+
/*
// gnuplot output for vitesse riemann
std::ofstream fout("ur_no_split");
@@ -462,6 +468,7 @@ public:
fout << xr[j][0] << ' ' << ur[j][0] << '\n';
}
*/
+
}
};