diff --git a/src/main.cpp b/src/main.cpp
index 8a6688090974596fc73e66e0dd69aa3d27114df1..2ab1f8d09538c6c6b87086930b462f9cfa57666b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -144,7 +144,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.8;
+ const double tmax=1.5;
double t=0.;
int itermax=std::numeric_limits<int>::max();
@@ -282,7 +282,7 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
// NAVIER-STOKES AVEC SPLITTING
-
+ /*
// Etape 1 du splitting - Euler
double dt_euler = 0.5*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
@@ -308,7 +308,7 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
-
+ */
// Diffusion pure
/*
double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
@@ -320,14 +320,14 @@ int main(int argc, char *argv[])
*/
// NAVIER-STOKES SANS SPLITTING
- /*
- double dt = 0.1*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+
+ double dt = 0.5*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
if (t+dt > tmax) {
dt = tmax-t;
}
no_splitting.computeNextStep(t,dt, unknowns);
t += dt;
- */
+
block_eos.updatePandCFromRhoE();
@@ -659,11 +659,11 @@ int main(int argc, char *argv[])
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");
+ std::ofstream fout("rho_3");
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';
}
}
@@ -691,16 +691,16 @@ int main(int argc, char *argv[])
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");
+ std::ofstream fout("u_3");
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(-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]*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';
}
}
@@ -715,10 +715,10 @@ int main(int argc, char *argv[])
//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] << ' ' << (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/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index a2bdcee9b41411e6208e07055cbc696d849527b3..a9219164d2dd3d736bac2de2595a326376288eae 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -73,7 +73,7 @@ public:
return m_xmax;
}
-
+ /*
// pas constant
Mesh(const Connectivity& connectivity)
@@ -82,8 +82,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();
@@ -91,7 +91,7 @@ public:
m_xr[r][0] = a + r*delta_x;
});
}
-
+ */
// pas non constant
/*
@@ -125,15 +125,15 @@ public:
// pas aleatoire
- /*
+
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes()),
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 h = (b-a)/connectivity.numberOfCells();
@@ -160,7 +160,7 @@ public:
//std::exit(0);
}
- */
+
~Mesh()
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 96278145bda4f445b14d73fbc2864671bccdb72a..60f4ddb6baa20d842c42df4665417cb7d3b9ef0b 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.014;
+ m_kj[j] = 0.5;
// TAC
@@ -392,11 +392,11 @@ 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.014*m_kj[j];
-
+ */
});
@@ -408,8 +408,8 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
- m_kL[0] = 0.;
- m_kR[0] = 0.;
+ m_kL[0] = 0.5;
+ m_kR[0] = 0.5;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index 6e97a84a1820819cdaabbdd34e17e54e6d4696ab..c1a9d64710c1c4d7a245ca0544c8aeefac01b0a5 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -190,19 +190,19 @@ private:
m_ur[r]=invAr(r)*br(r);
});
-
+ /*
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;
}
@@ -348,7 +348,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)));
@@ -374,7 +374,8 @@ public:
rhoj[j] = mj[j]/Vj[j];
});
- // Calcul de PT
+ /*
+ // Calcul de PT (1er essai)
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
double sum = 0;
double sum1 = 0;
@@ -382,24 +383,56 @@ public:
sum += (uj(cell_nodes(j,m)), Cjr(cell_nodes(j,m), m));
sum1 += Vj(cell_nodes(j,m));
}
-
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)*2.*sum/sum1;
}
+ */
+
+ // Calcul de PT (2eme essai, symetrisation)
+ const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
+ const Kokkos::View<const unsigned short*> face_nb_cells
+ = m_connectivity.faceNbCells();
+ const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces();
+ const Kokkos::View<const unsigned short*> cell_nb_faces
+ = m_connectivity.cellNbFaces();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+
+ std::vector<double> stock(2);
+ for (int l=0; l<cell_nb_faces(j); ++l) {
+ double sum = 0;
+ double sum2 = 0;
+ int k = cell_faces(j,l);
+ for (int i=0; i<face_nb_cells(k); ++i) {
+ int cell_here = face_cells(k,i);
+ sum += (1./Vj(cell_here))*uj[cell_here][0];
+ sum2 += 1./Vj(cell_here);
+ }
+ 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));
+ } 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);
+ } else {
+ //PTj(j) = pj(j) - kj(j)*(stock[1]-stock[0])/Vj(j);
+ PTj(j) = pj(j) - kj(j)*(stock[1]-stock[0])/Vj(j);
+ }
+
});
-
+
// Mise a jour de k
/*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0];
-
});
*/
}