plus de SIGFPE, mais solution numerique mauvaise

src/main.cpp

@@ -320,7 +320,7 @@ int main(int argc, char *argv[])
 
       // NAVIER-STOKES SANS SPLITTING
 
-      double dt = 0.05*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+      double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
       if (t+dt > tmax) {
 	dt = tmax-t;
       }
@@ -661,8 +661,8 @@ int main(int argc, char *argv[])
      std::ofstream fout("rho");
      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';
      }
      }
     
@@ -696,10 +696,10 @@ int main(int argc, char *argv[])
 
        //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';
      }
      }
 
@@ -714,10 +714,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';
      }
      }
 

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -352,8 +352,8 @@ 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.;
+  	//m_kj[j] =  xj[j][0];
+	m_kj[j] = 0.5;
 
 	// TAC
 	
@@ -408,8 +408,8 @@ public:
 
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 0.;
-    m_kR[0] = 1.;
+    m_kL[0] = 0.5;
+    m_kR[0] = 0.5;
 
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){  	

src/scheme/NoSplitting.hpp

@@ -312,6 +312,8 @@ public:
     Kokkos::View<double*> kj = unknowns.kj();
     Kokkos::View<double*> nuj = unknowns.nuj();
     Kokkos::View<double*> PTj = unknowns.PTj();
+    Kokkos::View<double*> kL = unknowns.kL();
+    Kokkos::View<double*> kR = unknowns.kR();
 
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
@@ -327,9 +329,17 @@ public:
 	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));
+	  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
@@ -350,6 +360,10 @@ public:
 	}
 	uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
 	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)));
+
       });
 
     // Calcul de e par la formule e = E-0.5 u^2
@@ -373,12 +387,12 @@ public:
       });
 
     // Mise a jour de k
-    
+    /*
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	kj(j) = xj[j][0];
 
       });
-    
+    */
   }
 };