From 6dd246ebe5191830b255c2fb55d1a80172b2198c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 26 Mar 2018 23:12:23 +0200
Subject: [PATCH] Fixed acoustic solver
---
experimental/AcousticSolver.cpp | 67 ++++++++++++++++++---------------
experimental/AcousticSolver.hpp | 6 ++-
main.cpp | 14 +++----
3 files changed, 47 insertions(+), 40 deletions(-)
diff --git a/experimental/AcousticSolver.cpp b/experimental/AcousticSolver.cpp
index f7ffd99b2..fc875ffad 100644
--- a/experimental/AcousticSolver.cpp
+++ b/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,28 +259,32 @@ 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]);
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+ 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;
Ej=new_Ej;
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
});
diff --git a/experimental/AcousticSolver.hpp b/experimental/AcousticSolver.hpp
index c8dde69ff..ba1462770 100644
--- a/experimental/AcousticSolver.hpp
+++ b/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:
diff --git a/main.cpp b/main.cpp
index c664fab2b..981358f6c 100644
--- a/main.cpp
+++ b/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();
--
GitLab