corrections de E et G pour point fixe ur

src/main.cpp

@@ -148,7 +148,7 @@ int main(int argc, char *argv[])
     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;
       }

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)
@@ -342,7 +345,6 @@ 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'; 
@@ -481,6 +500,18 @@ public:
       });
     */
     
+    // 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'; 
+    } 
+    */
   }
 };