diff --git a/experimental/AcousticSolver.cpp b/experimental/AcousticSolver.cpp index fc875ffadb5a48fc9150c28150e7eeb723d2d00e..16bf845f3233a2cc163a6183b123d7afeb0ce9b6 100644 --- a/experimental/AcousticSolver.cpp +++ b/experimental/AcousticSolver.cpp @@ -42,62 +42,43 @@ public: dst= Kokkos::reduction_identity<value_type>::min(); } }; - KOKKOS_INLINE_FUNCTION -double AcousticSolver:: -acoustic_dt(const Kokkos::View<const double*>& Vj, - const Kokkos::View<const double*>& cj) const +const Kokkos::View<const double*> +AcousticSolver::computeRhoCj(const Kokkos::View<const double*>& rhoj, + const Kokkos::View<const double*>& cj) { - const size_t nj = Vj.size(); - double dt = std::numeric_limits<double>::max(); - - Kokkos::View<double*> Vj_cj("Vj_cj", nj); - - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - Vj_cj[j] = Vj[j]/cj[j]; + Kokkos::View<double*> rhocj("rho_c", m_nj); + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + rhocj[j] = rhoj[j]*cj[j]; }); - - Kokkos::parallel_reduce(nj, ReduceMin(Vj_cj), dt); - - return dt; + return rhocj; } - KOKKOS_INLINE_FUNCTION -void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr, - const Kokkos::View<const double*>& xj, - const Kokkos::View<const double*>& rhoj, - const Kokkos::View<const double*>& uj, - 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, - Kokkos::View<double*>& ur, - Kokkos::View<double*[2]>& Fjr) const +const Kokkos::View<const double*[2]> +AcousticSolver::computeAjr(const Kokkos::View<const double*>& rhocj, + const Kokkos::View<const double*[2]>& Cjr) { - // calcul de ur - const size_t nr = ur.size(); - const size_t nj = uj.size(); - - Kokkos::View<double*> rhocj("rho_c", nj); - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { - rhocj[j] = rhoj[j]*cj[j]; - }); - - Kokkos::View<double*[2]> Ajr("Ajr", nj); - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { + Kokkos::View<double*[2]> Ajr("Ajr", m_nj); + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { for (int r=0; r<2; ++r) { Ajr(j,r) = rhocj(j)*Cjr(j,r)*Cjr(j,r); } }); - Kokkos::View<double*> Ar("Ar", nr); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) { + return Ajr; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolver::computeAr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const int*[2]>& node_cells, + const Kokkos::View<const int*[2]>& node_cell_local_node, + const Kokkos::View<const int*>& node_nb_cells) +{ + Kokkos::View<double*> Ar("Ar", m_nr); + Kokkos::parallel_for(m_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); @@ -107,14 +88,21 @@ 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); - }); - + return Ar; +} - Kokkos::View<double*> br("br", nr); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) { +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolver::computeBr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const double*[2]>& Cjr, + const Kokkos::View<const double*>& uj, + const Kokkos::View<const double*>& pj, + const Kokkos::View<const int*[2]>& node_cells, + const Kokkos::View<const int*[2]>& node_cell_local_node, + const Kokkos::View<const int*>& node_nb_cells) +{ + Kokkos::View<double*> br("br", m_nr); + Kokkos::parallel_for(m_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); @@ -124,69 +112,158 @@ void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr br(r) = sum; }); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) { + return br; +} + +KOKKOS_INLINE_FUNCTION +Kokkos::View<double*> +AcousticSolver::computeUr(const Kokkos::View<const double*>& Ar, + const Kokkos::View<const double*>& br) +{ + const Kokkos::View<const double*> invAr = inverse(Ar); + Kokkos::View<double*> ur("ur", m_nr); + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) { ur[r]=br(r)*invAr(r); }); ur[0]=0; - ur[nr-1]=0; + ur[m_nr-1]=0; - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { + return ur; +} + +KOKKOS_INLINE_FUNCTION +Kokkos::View<double*[2]> +AcousticSolver::computeFjr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const double*>& ur, + const Kokkos::View<const double*[2]>& Cjr, + const Kokkos::View<const double*>& uj, + const Kokkos::View<const double*>& pj, + const Kokkos::View<const int*[2]>& cell_nodes) +{ + Kokkos::View<double*[2]> Fjr("Fjr", m_nj); + + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { 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); } }); + + return Fjr; } -AcousticSolver::AcousticSolver(const long int& nj) +KOKKOS_INLINE_FUNCTION +double AcousticSolver:: +acoustic_dt(const Kokkos::View<const double*>& Vj, + const Kokkos::View<const double*>& cj) const { - Kokkos::View<double*> xj("xj", nj); - Kokkos::View<double*> rhoj("rhoj", nj); + double dt = std::numeric_limits<double>::max(); - Kokkos::View<double*> uj("uj", nj); + Kokkos::View<double*> Vj_cj("Vj_cj", m_nj); - Kokkos::View<double*> Ej("Ej", nj); - Kokkos::View<double*> ej("ej", nj); - Kokkos::View<double*> pj("pj", nj); - Kokkos::View<double*> Vj("Vj", nj); - Kokkos::View<double*> gammaj("gammaj", nj); - Kokkos::View<double*> cj("cj", nj); - Kokkos::View<double*> mj("mj", nj); - Kokkos::View<double*> inv_mj("inv_mj", nj); + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ + Vj_cj[j] = Vj[j]/cj[j]; + }); + + Kokkos::parallel_reduce(m_nj, ReduceMin(Vj_cj), dt); - const int nr=nj+1; + return dt; +} - Kokkos::View<double*> xr("xr", nr); +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolver::inverse(const Kokkos::View<const double*>& x) const +{ + Kokkos::View<double*> inv_x("inv_x", m_nr); + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) { + inv_x(r) = 1./x(r); + }); - const double delta_x = 1./nj; + return inv_x; +} - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ +KOKKOS_INLINE_FUNCTION +void AcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr, + const Kokkos::View<const double*>& xj, + const Kokkos::View<const double*>& rhoj, + const Kokkos::View<const double*>& uj, + 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, + Kokkos::View<double*>& ur, + Kokkos::View<double*[2]>& Fjr) +{ + const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj); + const Kokkos::View<const double*[2]> Ajr = computeAjr(rhocj, Cjr); + + const Kokkos::View<const double*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells); + const Kokkos::View<const double*> br = computeBr(Ajr, Cjr, uj, pj, + node_cells, node_cell_local_node, node_nb_cells); + + ur = computeUr(Ar, br); +#warning understand why this call diminishes performances + // Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes); + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + 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); + } + }); +} + +AcousticSolver::AcousticSolver(const long int& nj) + : m_nj(nj), + m_nr(nj+1) +{ + Kokkos::View<double*> xj("xj",m_nj); + Kokkos::View<double*> rhoj("rhoj",m_nj); + + Kokkos::View<double*> uj("uj",m_nj); + + Kokkos::View<double*> Ej("Ej",m_nj); + Kokkos::View<double*> ej("ej",m_nj); + Kokkos::View<double*> pj("pj",m_nj); + Kokkos::View<double*> Vj("Vj",m_nj); + Kokkos::View<double*> gammaj("gammaj",m_nj); + Kokkos::View<double*> cj("cj",m_nj); + Kokkos::View<double*> mj("mj",m_nj); + Kokkos::View<double*> inv_mj("inv_mj",m_nj); + + Kokkos::View<double*> xr("xr", m_nr); + + const double delta_x = 1./m_nj; + + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ xr[r] = r*delta_x; }); - Kokkos::View<int*[2]> cell_nodes("cell_nodes",nj,2); - Kokkos::View<int*[2]> node_cells("node_cells",nr,2); - Kokkos::View<int*[2]> node_cell_local_node("node_cells",nr,2); - Kokkos::View<int*> node_nb_cells("node_cells",nr); + Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2); + Kokkos::View<int*[2]> node_cells("node_cells",m_nr,2); + Kokkos::View<int*[2]> node_cell_local_node("node_cells",m_nr,2); + Kokkos::View<int*> node_nb_cells("node_cells",m_nr); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ node_nb_cells(r) = 2; }); node_nb_cells(0) = 1; - node_nb_cells(nr-1) = 1; + node_nb_cells(m_nr-1) = 1; node_cells(0,0) = 0; - Kokkos::parallel_for(nr-2, KOKKOS_LAMBDA(const int& r){ + Kokkos::parallel_for(m_nr-2, KOKKOS_LAMBDA(const int& r){ node_cells(r+1,0) = r; node_cells(r+1,1) = r+1; }); - node_cells(nr-1,0) = nj-1; + node_cells(m_nr-1,0) =m_nj-1; Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ cell_nodes(j,0) = j; cell_nodes(j,1) = j+1; }); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ for (int J=0; J<node_nb_cells(r); ++J) { int j = node_cells(r,J); for (int R=0; R<2; ++R) { @@ -243,7 +320,7 @@ AcousticSolver::AcousticSolver(const long int& nj) int itermax=std::numeric_limits<int>::max(); int iteration=0; - Kokkos::View<double*[2]> Cjr("Cjr", nj); + Kokkos::View<double*[2]> Cjr("Cjr",m_nj); Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { Cjr(j,0)=-1; Cjr(j,1)= 1; @@ -258,16 +335,17 @@ AcousticSolver::AcousticSolver(const long int& nj) t=tmax; } - Kokkos::View<double*> ur("ur", nr); - Kokkos::View<double*[2]> Fjr("Fjr", nr); - + Kokkos::View<double*> ur; + // Kokkos::View<double*[2]> Fjr; + Kokkos::View<double*[2]> Fjr("Fjr", m_nj); + computeExplicitFluxes(xr, xj, 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::View<double*> new_uj("new_uj",m_nj); + // Kokkos::View<double*> new_Ej("new_Ej",m_nj); Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ double momentum_fluxes = 0; @@ -277,25 +355,23 @@ AcousticSolver::AcousticSolver(const long int& nj) 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[j] -= dt/mj[j]*(momentum_fluxes); + Ej[j] -= dt/mj[j]*(energy_fluxes); }); - uj=new_uj; - Ej=new_Ej; + // uj=new_uj; + // Ej=new_Ej; Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { ej[j] = Ej[j] - 0.5 * uj[j]*uj[j]; }); - Kokkos::View<double*> xr_new("new_xr", nr); + // Kokkos::View<double*> xr_new("new_xr", m_nr); - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ - xr_new[r] = xr[r] + dt*ur[r]; + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ + xr[r] += dt*ur[r]; }); - xr = xr_new; - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ xj[j] = 0.5*(xr[j]+xr[j+1]); Vj[j] = xr[j+1]-xr[j]; @@ -310,6 +386,13 @@ AcousticSolver::AcousticSolver(const long int& nj) ++iteration; } + { + std::ofstream fout("rho"); + for (int j=0; j<nj; ++j) { + fout << xj[j] << ' ' << rhoj[j] << '\n'; + } + } + std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; diff --git a/experimental/AcousticSolver.hpp b/experimental/AcousticSolver.hpp index ba1462770b4bb1018f0f4945cfe0892bce2a09ab..d8429d5acc023acea9ef6930c8de1306f0bc1370 100644 --- a/experimental/AcousticSolver.hpp +++ b/experimental/AcousticSolver.hpp @@ -7,8 +7,43 @@ class AcousticSolver { private: - inline double e(double rho, double p, double gamma) const; - inline double p(double rho, double e, double gamma) const; + inline const Kokkos::View<const double*> + computeRhoCj(const Kokkos::View<const double*>& rhoj, + const Kokkos::View<const double*>& cj); + + inline const Kokkos::View<const double*[2]> + computeAjr(const Kokkos::View<const double*>& rhocj, + const Kokkos::View<const double*[2]>& Cjr); + + inline const Kokkos::View<const double*> + computeAr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const int*[2]>& node_cells, + const Kokkos::View<const int*[2]>& node_cell_local_node, + const Kokkos::View<const int*>& node_nb_cells); + + inline const Kokkos::View<const double*> + computeBr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const double*[2]>& Cjr, + const Kokkos::View<const double*>& uj, + const Kokkos::View<const double*>& pj, + const Kokkos::View<const int*[2]>& node_cells, + const Kokkos::View<const int*[2]>& node_cell_local_node, + const Kokkos::View<const int*>& node_nb_cells); + + Kokkos::View<double*> + computeUr(const Kokkos::View<const double*>& Ar, + const Kokkos::View<const double*>& br); + + Kokkos::View<double*[2]> + computeFjr(const Kokkos::View<const double*[2]>& Ajr, + const Kokkos::View<const double*>& ur, + const Kokkos::View<const double*[2]>& Cjr, + const Kokkos::View<const double*>& uj, + const Kokkos::View<const double*>& pj, + const Kokkos::View<const int*[2]>& cell_nodes); + + const Kokkos::View<const double*> + inverse(const Kokkos::View<const double*>& x) const; inline double acoustic_dt(const Kokkos::View<const double*>& Vj, const Kokkos::View<const double*>& cj) const; @@ -26,9 +61,13 @@ private: const Kokkos::View<const int*>& node_nb_cells, const Kokkos::View<const int*[2]>& node_cell_local_node, Kokkos::View<double*>& ur, - Kokkos::View<double*[2]>& Fjr) const; + Kokkos::View<double*[2]>& Fjr); struct ReduceMin; + + const long int m_nj; + const long int m_nr; + public: AcousticSolver(const long int& nj); };