diff --git a/src/main.cpp b/src/main.cpp
index 19857301af9b35b4d57b7509795b9eb75e8a57dc..e16fdcc70f84a7b2903b4f19ad52f688fe2f3c60 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -148,7 +148,6 @@ int main(int argc, char *argv[])
double t=0.;
int itermax=std::numeric_limits<int>::max();
- //int itermax = 202;
int iteration=0;
Kokkos::View<double*> rhoj = unknowns.rhoj();
@@ -302,7 +301,7 @@ int main(int argc, char *argv[])
// NAVIER-STOKES AVEC SPLITTING
/*
// Etape 1 du splitting - Euler
- double dt_euler = 0.1*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
}
@@ -310,7 +309,7 @@ int main(int argc, char *argv[])
t += dt_euler;
// Etape 2 du splitting - Diffusion
- double param = 1.;
+ double param = 0.9;
double dt_diff = param*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
double t_diff = t-dt_euler;
if (dt_euler <= dt_diff) {
@@ -350,6 +349,7 @@ int main(int argc, char *argv[])
}
no_splitting.computeNextStep(t,dt, unknowns);
t += dt;
+
/*
// Fichier pas de temps
std::ofstream past("pas_no_split", std::ios::app);
@@ -362,7 +362,6 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
- //std::exit(0);
// Animations
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 03cf62de92d8da737bdd4a45f5fec418a09cc891..40298dc5c01217cc5002c3f14f78e12559406ffb 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -181,19 +181,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;
}
@@ -379,13 +379,6 @@ public:
nuj(j) = 0.5*(1.+xj[j][0]);
});
*/
-
- // gnuplot output for vitesse riemann
- std::ofstream fout("ur_split");
- fout.precision(15);
- for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
- fout << xr[j][0] << ' ' << ur[j][0] << '\n';
- }
}
};
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 281000de92e68e564f57c1ce8b7604c78bb4d406..dd3b42900f4edde1c4d98c33a520cf103f782dbd 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)));
@@ -123,13 +123,13 @@ 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) = -((kR(0)/Vj(m_mesh.numberOfCells()-1) + kj(cell_here)/Vj(cell_here))/(1./Vj(cell_here) + 1./Vj(m_mesh.numberOfCells()-1)))*(1./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
// k = 0.5
- m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
- m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
+ //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
+ //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.));
@@ -179,16 +179,16 @@ private:
});
// Conditions aux bords
- /*
+
m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
- */
- // Kidder
+ // 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 ;
@@ -495,7 +495,7 @@ public:
//Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + xj[j][0]*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]) + (1.+xj[j][0])*((3.*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0])-xj[j][0]*pi*pi*pi*pi*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])))*(std::exp(-2.*t)-1.)+0.5*pi*pi*std::exp(-2.*t)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
// ajout second membre pour kidder (k = 0.5)
- Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ //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));
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index c441063fe62f73c167cadc7c2a216b5d51391f85..f25ddc0a631db2c9347c92698d27e118749b8122 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -392,53 +392,45 @@ public:
fout3 << xj[j][0] << ' ' << PTj[j] << '\n';
}
*/
+
// Calcul de PT (3eme essai, avec uR du solveur de Riemann)
- 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 {
+ // iteration pour point fixe ur
+ for (int itconv=0; itconv<100; ++itconv){
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
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);
- // }
- });
- // 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) {
- // });
-
+ });
+
+ // Copie
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur0[inode][0]=m_ur[inode][0];
+ }
- // 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;
+ // Calcule les flux
+ computeExplicitFluxes(xr, xj, rhoj, kj, uj, PTj, cj, Vj, Cjr, t);
+
+ // Relaxation de ur
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur[inode][0]=0.5*m_ur[inode][0]+0.5*m_ur0[inode][0];
+ }
+
+ // Calcul erreur L1
+ 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) {
@@ -478,7 +470,7 @@ 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);
@@ -487,31 +479,19 @@ public:
}
// gnuplot output for vitesse riemann
- std::ofstream fout1("ur202");
+ 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];
});
*/
-
- // 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';
- }
- */
+
}
};