Fixed acoustic solver

experimental/AcousticSolver.cpp

@@ -72,11 +72,13 @@ void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr
 					   const Kokkos::View<const double*>& pj,
 					   const Kokkos::View<const double*>& cj,
 					   const Kokkos::View<const double*>& Vj,
+					   const Kokkos::View<const double*[2]>& Cjr,
+					   const Kokkos::View<const int*[2]>& cell_nodes,
 					   const Kokkos::View<const int*[2]>& node_cells,
 					   const Kokkos::View<const int*>& node_nb_cells,
-					   const Kokkos::View<const int*[2]> node_cell_local_node,
+					   const Kokkos::View<const int*[2]>& node_cell_local_node,
 					   Kokkos::View<double*>& ur,
-					   Kokkos::View<double*>& pr) const
+					   Kokkos::View<double*[2]>& Fjr) const
 {
   // calcul de ur
   const size_t nr = ur.size();
@@ -87,13 +89,6 @@ void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr
       rhocj[j] = rhoj[j]*cj[j];
     });
 
-
-  Kokkos::View<double*[2]> Cjr("Cjr", nj);
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-      Cjr(j,0)=-1;
-      Cjr(j,1)= 1;
-    });
-
   Kokkos::View<double*[2]> Ajr("Ajr", nj);
   Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
       for (int r=0; r<2; ++r) {
@@ -101,7 +96,6 @@ void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr
       }
     });
 
-
   Kokkos::View<double*> Ar("Ar", nr);
   Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) {
       double sum = 0;
@@ -113,32 +107,33 @@ void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr
       Ar(r) = sum;
     });
 
+  Kokkos::View<double*> invAr("invAr", nr);
+  Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) {
+      invAr(r) = 1./Ar(r);
+    });
+
+
   Kokkos::View<double*> br("br", nr);
   Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) {
       double sum = 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);
-  	sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj[j];
+  	sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj[J];
       }
       br(r) = sum;
     });
 
   Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) {
-      ur[r]=br(r)/Ar(r);
+      ur[r]=br(r)*invAr(r);
     });
   ur[0]=0;
   ur[nr-1]=0;
 
-  // calcul de pr
-  pr[0] = pj[0] + rhoj[0]*cj[0]*(ur[0] - uj[0]);
   Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
-    const int r = j+1;
-
-    const double ujr = uj[j];
-    const double pjr = pj[j];
-
-    pr[r]=pjr+rhoj[j]*cj[j]*(ujr-ur[r]);
+      for (int r=0; r<2; ++r) {
+	Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+Cjr(j,r)*pj(j);
+      }
     });
 }
 
@@ -248,6 +243,12 @@ AcousticSolver::AcousticSolver(const long int& nj)
   int itermax=std::numeric_limits<int>::max();
   int iteration=0;
 
+  Kokkos::View<double*[2]> Cjr("Cjr", nj);
+  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
+      Cjr(j,0)=-1;
+      Cjr(j,1)= 1;
+    });
+
   while((t<tmax) and (iteration<itermax)) {
     double dt = 0.4*acoustic_dt(Vj, cj);
     if (t+dt<tmax) {
@@ -258,22 +259,26 @@ AcousticSolver::AcousticSolver(const long int& nj)
     }
     
     Kokkos::View<double*> ur("ur", nr);
-    Kokkos::View<double*> pr("pr", nr);
+    Kokkos::View<double*[2]> Fjr("Fjr", nr);
 
     computeExplicitFluxes(xr, xj,
-			  rhoj, uj, pj, cj, Vj,
-			  node_cells, node_nb_cells, node_cell_local_node,
-			  ur, pr);
+			  rhoj, uj, pj, cj, Vj, Cjr,
+			  cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
+			  ur, Fjr);
 
     Kokkos::View<double*> new_uj("new_uj", nj);
     Kokkos::View<double*> new_Ej("new_Ej", nj);
 
     Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-	int rm=cell_nodes(j,0);
-	int rp=cell_nodes(j,1);
-
-	new_uj[j] = uj[j] + dt/mj[j]*(pr[rm]-pr[rp]);
-	new_Ej[j] = Ej[j] + dt/mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]);
+	double momentum_fluxes = 0;
+	double energy_fluxes = 0;
+	for (int R=0; R<2; ++R) {
+	  const int r=cell_nodes(j,R);
+	  momentum_fluxes += Fjr(j,R);
+	  energy_fluxes += Fjr(j,R)*ur[r];
+	}
+	new_uj[j] = uj[j] - dt/mj[j]*(momentum_fluxes);
+	new_Ej[j] = Ej[j] - dt/mj[j]*(energy_fluxes);
       });
 
     uj=new_uj;

experimental/AcousticSolver.hpp

@@ -20,11 +20,13 @@ private:
 				    const Kokkos::View<const double*>& pj,
 				    const Kokkos::View<const double*>& cj,
 				    const Kokkos::View<const double*>& Vj,
+				    const Kokkos::View<const double*[2]>& Cjr,
+				    const Kokkos::View<const int*[2]>& cell_nodes,
 				    const Kokkos::View<const int*[2]>& node_cells,
 				    const Kokkos::View<const int*>& node_nb_cells,
-				    const Kokkos::View<const int*[2]> node_cell_local_node,
+				    const Kokkos::View<const int*[2]>& node_cell_local_node,
 				    Kokkos::View<double*>& ur,
-				    Kokkos::View<double*>& pr) const;
+				    Kokkos::View<double*[2]>& Fjr) const;
 
   struct ReduceMin;
 public:

main.cpp

@@ -91,13 +91,13 @@ int main(int argc, char *argv[])
     MeshLessAcousticSolver acoustic_solver(number);
     method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
   }
-#warning UNCOMMENT TO CONTINUE
-  // { // class for acoustic solver
-  //   Kokkos::Timer timer;
-  //   timer.reset();
-  //   AcousticSolver acoustic_solver(number);
-  //   method_cost_map["AcousticSolver"] = timer.seconds();
-  // }
+
+  { // class for acoustic solver
+    Kokkos::Timer timer;
+    timer.reset();
+    AcousticSolver acoustic_solver(number);
+    method_cost_map["AcousticSolver"] = timer.seconds();
+  }
 
   Kokkos::finalize();