diff --git a/src/main.cpp b/src/main.cpp
index e2781da1d09b26ebb2d7b51194036b5841b65ab5..19857301af9b35b4d57b7509795b9eb75e8a57dc 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -143,12 +143,12 @@ 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;
double t=0.;
int itermax=std::numeric_limits<int>::max();
- //int itermax = 500;
+ //int itermax = 202;
int iteration=0;
Kokkos::View<double*> rhoj = unknowns.rhoj();
@@ -341,7 +341,7 @@ int main(int argc, char *argv[])
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);
- double dt_diff = 0.1*no_splitting.nosplitting_dt(Vj,cj,rhoj,kj);
+ double dt_diff = 0.9*no_splitting.nosplitting_dt(Vj,cj,rhoj,kj);
if (dt_diff < dt) {
dt = dt_diff;
}
@@ -362,7 +362,7 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
- // std::exit(0);
+ //std::exit(0);
// Animations
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index 25345194e8909ed656c1a0f3ef311b72dc570ef8..c441063fe62f73c167cadc7c2a216b5d51391f85 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -248,6 +248,7 @@ private:
Kokkos::View<Rdd*> m_inv_Ar;
Kokkos::View<Rd**> m_Fjr;
Kokkos::View<Rd*> m_ur;
+ Kokkos::View<Rd*> m_ur0;
Kokkos::View<double*> m_rhocj;
Kokkos::View<double*> m_Vj_over_cj;
@@ -263,6 +264,7 @@ public:
m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
m_ur("ur", m_mesh.numberOfNodes()),
+ m_ur0("ur0", m_mesh.numberOfNodes()),
m_rhocj("rho_c", m_mesh.numberOfCells()),
m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
{
@@ -320,6 +322,7 @@ public:
const Kokkos::View<const Rd**> Fjr = m_Fjr;
const Kokkos::View<const Rd*> ur = m_ur;
+ const Kokkos::View<const Rd*> ur0 = m_ur0;
/*
// Calcul de PT (1er essai)
@@ -341,8 +344,7 @@ public:
}
});
- */
-
+ */
/*
// Calcul de PT (2eme essai, symetrisation)
const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
@@ -377,7 +379,7 @@ public:
}
});
- */
+
std::ofstream fout2("pj");
fout2.precision(15);
@@ -389,37 +391,54 @@ public:
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, kj, uj, pj, cj, Vj, Cjr, t);
+ for (int itconv=0; itconv<100; ++itconv){
+
+ // computeExplicitFluxes(xr, xj, rhoj, kj, 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)*(ur[cell_nodes(j,1)][0])/Vj(j);
- } else if (j == m_mesh.numberOfCells()-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 {
+ // 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)*(ur[cell_nodes(j,1)][0])/Vj(j);
+ // } else if (j == m_mesh.numberOfCells()-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);
sum += (ur(node_here), Cjr(j,k));
}
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);
- }
+ // }
});
+ // jespere que ca copie, ca...
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur0[inode][0]=m_ur[inode][0];
+ }
+ // m_ur0=m_ur;
+ // Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ // });
+
-
// Calcule les flux
computeExplicitFluxes(xr, xj, rhoj, kj, uj, PTj, cj, Vj, Cjr, t);
-
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur[inode][0]=0.7*m_ur[inode][0]+0.3*m_ur0[inode][0];
+ }
+ double sum=0.;
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ sum+=std::abs(m_ur0[inode][0]-m_ur[inode][0]);
+ }
+ sum/=double(m_mesh.numberOfNodes());
+ std::cout << " it " << itconv << " sum " << sum << std::endl;
+ if(sum<1.e-6) break;
+ }
// 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) {
@@ -468,7 +487,7 @@ public:
}
// gnuplot output for vitesse riemann
- std::ofstream fout1("ur");
+ std::ofstream fout1("ur202");
fout1.precision(15);
for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
@@ -480,7 +499,19 @@ public:
kj(j) = xj[j][0];
});
*/
-
+
+ // stocke la vitesse pour la prochaine iteration
+ // Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
+ // ur0[r][0] = ur[r][0];
+ // });
+ /*
+ // gnuplot output for vitesse riemann
+ std::ofstream fout1("ur0");
+ fou.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fou << xr[j][0] << ' ' << ur0[j][0] << '\n';
+ }
+ */
}
};