kidder convergence ok pour maillage unif mais pb avec maillage aleatoire

src/main.cpp

@@ -284,7 +284,7 @@ int main(int argc, char *argv[])
       // NAVIER-STOKES AVEC SPLITTING
       /*
       // Etape 1 du splitting - Euler
-      double dt_euler = 0.9*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;
       }

src/mesh/Mesh.hpp

@@ -75,7 +75,7 @@ public:
 
 
   // pas constant
-
+  /*
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes()),
@@ -91,7 +91,7 @@ public:
   	m_xr[r][0] = a + r*delta_x;
       });
   }
-  
+  */
 
   // pas non constant
   /*
@@ -125,7 +125,7 @@ public:
 
   
   // pas aleatoire
-  /*
+  
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes()),
@@ -143,21 +143,24 @@ public:
     std::mt19937 mt(rd());
     std::uniform_real_distribution<double> dist(-h/2.1,h/2.1); 
 
-    // creation fichier pour tracer h en fonction de x 
-    std::ofstream fout("alea");
-    
     for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
       const double delta_xr = dist(mt);
       m_xr[r][0] = r*h + delta_xr;
-      fout << m_xr[r][0] << ' ' << m_xr[r][0]-m_xr[r-1][0] << '\n';
-      //std::cout << m_xr[r][0] << std::endl;
     } 
    
-    fout.close();
+    // creation fichier pour tracer h en fonction de x 
+    //std::ofstream fout("alea");
+    //for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
+    //const double delta_xr = dist(mt);
+    //m_xr[r][0] = r*h + delta_xr;
+    //fout << m_xr[r][0] << ' ' << m_xr[r][0]-m_xr[r-1][0] << '\n';
+    //std::cout << m_xr[r][0] << std::endl;
+    //} 
+    //fout.close();
+    //std::exit(0);
    
-    std::exit(0);
   }
-  */
+
   
 
   ~Mesh()

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -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.005;
 
 	// TAC
 	
@@ -401,15 +401,15 @@ public:
       });
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_PTj(j) = 0;
+	m_PTj(j) = m_pj(j);
       });
     
      // Conditions aux bords de Dirichlet sur u et k
 
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 0.5;
-    m_kR[0] = 0.5;
+    m_kL[0] = 0.005;
+    m_kR[0] = 0.005;
 
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){  	

src/scheme/NoSplitting.hpp

@@ -314,9 +314,12 @@ public:
     Kokkos::View<double*> PTj = unknowns.PTj();
     Kokkos::View<double*> kL = unknowns.kL();
     Kokkos::View<double*> kR = unknowns.kR();
+    Kokkos::View<Rd*> uL = unknowns.uL();
+    Kokkos::View<Rd*> uR = unknowns.uR();
 
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+    const Kokkos::View<const double*> Vl = m_mesh_data.Vl();
     const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
     Kokkos::View<Rd*> xr = m_mesh.xr();
 
@@ -325,23 +328,6 @@ public:
     const Kokkos::View<const unsigned short*> cell_nb_nodes
       = m_connectivity.cellNbNodes();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-	double sum = 0;
-	for (int m=0; m<cell_nb_nodes(j); ++m) {
-	  //sum += uj[cell_nodes(j,m)][0];
-	  sum += (uj(cell_nodes(j,m)), Cjr(cell_nodes(j,m), m));
-	}
-
-	if (j == 0) {
-	  PTj(j) = pj(j) - kL(0)*sum/Vj(j);
-	} else if (j == m_mesh.numberOfCells()-1) {
-	  PTj(j) = pj(j) - kR(0)*sum/Vj(j);
-	} else {
-	  PTj(j) = pj(j) - kj(j)*sum/Vj(j);
-	}
-
-      });
-
     // Calcule les flux
     computeExplicitFluxes(xr, xj, kj, rhoj, uj, PTj, cj, Vj, Cjr, t);
 
@@ -363,7 +349,9 @@ public:
 
 	// 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)));
-
+	// 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))); 
       });
 
     // Calcul de e par la formule e = E-0.5 u^2
@@ -386,6 +374,28 @@ public:
 	rhoj[j] = mj[j]/Vj[j];
       });
 
+    // Calcul de PT
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	double sum = 0;
+	double sum1 = 0;
+	for (int m=0; m<cell_nb_nodes(j); ++m) {
+	  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));
+	} 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)*sum/Vl(j);
+	  PTj(j) = pj(j) - kj(j)*2.*sum/sum1;
+	}
+
+      });
+    
     // Mise a jour de k
     /*
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {