correction CFL et position R(t) a la place de x(t) dans Cl de computeUr

src/main.cpp

@@ -138,7 +138,7 @@ 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=1.;
+    const double tmax=0.7;
     double t=0;
 
     int itermax=std::numeric_limits<int>::max();
@@ -185,7 +185,6 @@ int main(int argc, char *argv[])
 	  }
 	  t_diff += dt_diff;
 	}
-	std::cout << "diff : " << t_diff << std::endl;
       }
       
 

src/mesh/Mesh.hpp

@@ -15,6 +15,8 @@ public:
 private:
   const Connectivity& m_connectivity;
   Kokkos::View<Rd*> m_xr;
+  Kokkos::View<Rd*> m_x0;
+  Kokkos::View<Rd*> m_xmax;
 
 public:
 
@@ -48,16 +50,43 @@ public:
     return m_xr;
   }
 
+  Kokkos::View<Rd*> x0()
+  {
+    return m_x0;
+  }
+
+  const Kokkos::View<const Rd*> x0() const
+  {
+    return m_x0;
+  }
+
+  Kokkos::View<Rd*> xmax()
+  {
+    return m_xmax;
+  }
+
+  const Kokkos::View<const Rd*> xmax() const
+  {
+    return m_xmax;
+  }
+
+
   // pas constant
 
   
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
-      m_xr("xr", connectivity.numberOfNodes())
+      m_xr("xr", connectivity.numberOfNodes()),
+      m_x0("x0", 1),
+      m_xmax("xmax", 1)
   {
-    const double delta_x = (1.-0.1)/connectivity.numberOfCells();
+    double a = 0.1;
+    double b = 1.;
+    m_x0[0][0] = a;
+    m_xmax[0][0] = b;
+    const double delta_x = (b-a)/connectivity.numberOfCells();
     Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
-  	m_xr[r][0] = 0.1+r*delta_x;
+  	m_xr[r][0] = a + r*delta_x;
       });
   }
   

src/scheme/AcousticSolver.hpp

@@ -129,10 +129,12 @@ private:
   computeBr(const Kokkos::View<const Rdd**>& Ajr,
 	    const Kokkos::View<const Rd**>& Cjr,
 	    const Kokkos::View<const Rd*>& uj,
-	    const Kokkos::View<const double*>& pj) {
+	    const Kokkos::View<const double*>& pj,
+	    const double t) {
     const Kokkos::View<const unsigned int**>& node_cells = m_connectivity.nodeCells();
     const Kokkos::View<const unsigned short**>& node_cell_local_node = m_connectivity.nodeCellLocalNode();
     const Kokkos::View<const unsigned short*>& node_nb_cells = m_connectivity.nodeNbCells();
+    Kokkos::View<Rd*> xr = m_mesh.xr();
 
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
 	Rd& br = m_br(r);
@@ -142,6 +144,22 @@ private:
 	  const int R = node_cell_local_node(r,j);
 	  br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
 	}
+	// if (r==0){
+	//   for (int j=0; j<node_nb_cells(r); ++j) {
+	//     const int J = node_cells(r,j);
+	//     const int R = node_cell_local_node(r,j); 
+	//     double rho0=3./(ht*ht)*(xr(0),xr(0)); 
+	//     double p0=2*rho0*rho0*rho0;
+	//     br-=p0*Cjr(J,R);
+	//   }
+	// }
+	// if (r==(m_mesh.numberOfNodes()-1)){
+	//   for (int j=0; j<node_nb_cells(r); ++j) {
+	//     const int J = node_cells(r,j);
+	//     const int R = node_cell_local_node(r,j);
+	//     br-=pmax*Cjr(J,R);
+	//   }
+	// }
       });
 
     return m_br;
@@ -172,14 +190,21 @@ private:
     const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
 
     Kokkos::View<Rd*> xr = m_mesh.xr();
+    Kokkos::View<Rd*> x0 = m_mesh.x0();
+    Kokkos::View<Rd*> xmax = m_mesh.xmax();
 
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
 	m_ur[r]=invAr(r)*br(r);
       });
 
     // Conditions aux limites dependant du temps
-    m_ur[0]=(-t/((50./9.)-t*t))*xr[0];
-    m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
+    //m_ur[0]=(-t/((50./9.)-t*t))*xr[0];
+    //m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
+
+    // R(t) = x*h(t) a la place de x(t) 
+    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;
   }
@@ -236,7 +261,7 @@ private:
     const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
 
     const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
-    const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj);
+    const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj, t);
 
     Kokkos::View<Rd*> ur = m_ur;
     Kokkos::View<Rd**> Fjr = m_Fjr;

src/scheme/FiniteVolumesDiffusion.hpp

@@ -124,14 +124,14 @@ private:
     //m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
     */
 
-    // k = x
-    m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
-    m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
-    
     // 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;
 
+    // k = x
+    m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
+    m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
+    
     return m_Fl ;
   }
 
@@ -264,7 +264,8 @@ public:
 	 sum += kj(cell_nodes(j,m));
 	}
 
-	dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
+	// dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
+	dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl;
 	
       });