diff --git a/experimental/AcousticSolver.cpp b/experimental/AcousticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f7ffd99b25f6b6ff4e262cb2bd4fcded808bc3de --- /dev/null +++ b/experimental/AcousticSolver.cpp @@ -0,0 +1,311 @@ +#include <AcousticSolver.hpp> +#include <rang.hpp> + +#include <memory> + +#include <BlockPerfectGas.hpp> + +typedef const double my_double; + +struct AcousticSolver::ReduceMin +{ +private: + const Kokkos::View<my_double*> x_; + +public: + typedef Kokkos::View<my_double*>::non_const_value_type value_type; + + ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {} + + typedef Kokkos::View<my_double*>::size_type size_type; + + KOKKOS_INLINE_FUNCTION void + operator() (const size_type i, value_type& update) const + { + if (x_(i) < update) { + update = x_(i); + } + } + + KOKKOS_INLINE_FUNCTION void + join (volatile value_type& dst, + const volatile value_type& src) const + { + if (src < dst) { + dst = src; + } + } + + KOKKOS_INLINE_FUNCTION void + init (value_type& dst) const + { // The identity under max is -Inf. + 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 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::parallel_reduce(nj, ReduceMin(Vj_cj), dt); + + return dt; +} + + +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 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*>& pr) const +{ + // 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]> 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) { + 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) { + 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); + } + Ar(r) = sum; + }); + + 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]; + } + br(r) = sum; + }); + + Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r) { + ur[r]=br(r)/Ar(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]); + }); +} + +AcousticSolver::AcousticSolver(const long int& nj) +{ + Kokkos::View<double*> xj("xj", nj); + Kokkos::View<double*> rhoj("rhoj", nj); + + Kokkos::View<double*> uj("uj", 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); + + const int nr=nj+1; + + Kokkos::View<double*> xr("xr", nr); + + const double delta_x = 1./nj; + + Kokkos::parallel_for(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::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ + node_nb_cells(r) = 2; + }); + node_nb_cells(0) = 1; + node_nb_cells(nr-1) = 1; + + node_cells(0,0) = 0; + Kokkos::parallel_for(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; + + 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){ + for (int J=0; J<node_nb_cells(r); ++J) { + int j = node_cells(r,J); + for (int R=0; R<2; ++R) { + if (cell_nodes(j,R) == r) { + node_cell_local_node(r,J) = R; + } + } + } + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + Vj[j] = xr[cell_nodes(j,1)]-xr[cell_nodes(j,0)]; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + if (xj[j]<0.5) { + rhoj[j]=1; + } else { + rhoj[j]=0.125; + } + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + if (xj[j]<0.5) { + pj[j]=1; + } else { + pj[j]=0.1; + } + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + gammaj[j] = 1.4; + }); + + BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); + + block_eos.updateEandCFromRhoP(); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + Ej[j] = ej[j]+0.5*uj[j]*uj[j]; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + mj[j] = rhoj[j] * Vj[j]; + }); + + const double tmax=0.2; + double t=0; + + int itermax=std::numeric_limits<int>::max(); + int iteration=0; + + while((t<tmax) and (iteration<itermax)) { + double dt = 0.4*acoustic_dt(Vj, cj); + if (t+dt<tmax) { + t+=dt; + } else { + dt=tmax-t; + t=tmax; + } + + Kokkos::View<double*> ur("ur", nr); + Kokkos::View<double*> pr("pr", nr); + + computeExplicitFluxes(xr, xj, + rhoj, uj, pj, cj, Vj, + node_cells, node_nb_cells, node_cell_local_node, + ur, pr); + + 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]); + }); + + 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::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ + xr_new[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]; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + rhoj[j] = mj[j]/Vj[j]; + }); + + block_eos.updatePandCFromRhoE(); + + ++iteration; + } + + 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 new file mode 100644 index 0000000000000000000000000000000000000000..c8dde69ff0dfeeaf434b2a9ed8d9585cb1818ecc --- /dev/null +++ b/experimental/AcousticSolver.hpp @@ -0,0 +1,34 @@ +#ifndef ACOUSTIC_SOLVER_HPP +#define ACOUSTIC_SOLVER_HPP + +#include <Kokkos_Core.hpp> + +class AcousticSolver +{ +private: + + inline double e(double rho, double p, double gamma) const; + inline double p(double rho, double e, double gamma) const; + + inline double acoustic_dt(const Kokkos::View<const double*>& Vj, + const Kokkos::View<const double*>& cj) const; + + inline void 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 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*>& pr) const; + + struct ReduceMin; +public: + AcousticSolver(const long int& nj); +}; + +#endif // ACOUSTIC_SOLVER_HPP diff --git a/experimental/CMakeLists.txt b/experimental/CMakeLists.txt index 52ff12b08290eceb152104d353591c003f3c7678..bd8990bd83911408a40c53a874945e753d732b33 100644 --- a/experimental/CMakeLists.txt +++ b/experimental/CMakeLists.txt @@ -8,6 +8,7 @@ include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include) add_library( PastisExperimental + AcousticSolver.cpp RawKokkosAcousticSolver.cpp MeshLessAcousticSolver.cpp ) diff --git a/main.cpp b/main.cpp index 7141464e0f5e3334a4dd6f53350fdaa992ee1585..c664fab2b14e088aaea21e8652e03fd08de0ad25 100644 --- a/main.cpp +++ b/main.cpp @@ -8,6 +8,7 @@ #include <RawKokkosAcousticSolver.hpp> #include <MeshLessAcousticSolver.hpp> +#include <AcousticSolver.hpp> #include <CLI/CLI.hpp> #include <cassert> @@ -76,19 +77,27 @@ int main(int argc, char *argv[]) Kokkos::DefaultExecutionSpace::print_configuration(std::cout); std::map<std::string, double> method_cost_map; - { + + { // Basic function based acoustic solver Kokkos::Timer timer; timer.reset(); RawKokkos::AcousticSolver(number); method_cost_map["RawKokkos"] = timer.seconds(); } - { + { // class for acoustic solver (mesh less) Kokkos::Timer timer; timer.reset(); 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(); + // } Kokkos::finalize(); @@ -98,7 +107,6 @@ int main(int argc, char *argv[]) for (const auto& method_cost : method_cost_map) { size = std::max(size, method_cost.first.size()); } - size+=2; for (const auto& method_cost : method_cost_map) { std::cout << "* [" << std::setw(size) << std::left