diff --git a/.gitignore b/.gitignore index d55e76cb535b79f68639000251f8833a9a291a26..234c9dd01b7bddde1ab85d0085fd6c2edb152e51 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +*.o +*.a *~ build/ CMakeFiles/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 7bcbe7fe9b1954684593b3202dfe365b369b9a3d..7c09d21d2eb538f2f7142e119b261c3b6518d3f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,24 +16,32 @@ project (Pastis set(PASTIS_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}") set(PASTIS_BINARY_DIR "${CMAKE_CURRENT_BINARY_DIR}") -# Rang (colors? Useless thus necessary!) -include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include) - -# CLI11 -include_directories(${PASTIS_SOURCE_DIR}/packages/CLI11/include) +#------------------------------------------------------ # Kokkos set(KOKKOS_ENABLE_OPENMP ON CACHE BOOL "") add_subdirectory(${PASTIS_SOURCE_DIR}/packages/kokkos) include_directories(${Kokkos_INCLUDE_DIRS_RET}) +# Compiler flags +include(GetKokkosCompilerFlags) +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall") + +#------------------------------------------------------ + +# Rang (colors? Useless thus necessary!) +include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include) + +# CLI11 +include_directories(${PASTIS_SOURCE_DIR}/packages/CLI11/include) + # Pastis utils add_subdirectory(utils) include_directories(utils) -# Compiler flags -include(GetKokkosCompilerFlags) -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall") +# Pastis experimental +add_subdirectory(experimental) +include_directories(experimental) # ---------------- Checks for includes ---------------- @@ -58,5 +66,6 @@ add_executable( target_link_libraries( pastis kokkos - PastisUtils) + PastisUtils + PastisExperimental) diff --git a/experimental/AcousticSolver.cpp b/experimental/AcousticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3088ff7176c03d46e6e429cb0b503123d0ecd66f --- /dev/null +++ b/experimental/AcousticSolver.cpp @@ -0,0 +1,385 @@ +#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 +const Kokkos::View<const double*> +AcousticSolver::computeRhoCj(const Kokkos::View<const double*>& rhoj, + const Kokkos::View<const double*>& cj) +{ + Kokkos::View<double*> rhocj("rho_c", m_nj); + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + rhocj[j] = rhoj[j]*cj[j]; + }); + return rhocj; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*[2]> +AcousticSolver::computeAjr(const Kokkos::View<const double*>& rhocj, + const Kokkos::View<const double*[2]>& Cjr) +{ + 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); + } + }); + + 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); + const int R = node_cell_local_node(r,j); + sum += Ajr(J,R); + } + Ar(r) = sum; + }); + + return Ar; +} + +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); + const int R = node_cell_local_node(r,j); + sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj[J]; + } + br(r) = sum; + }); + + 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[m_nr-1]=0; + + 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; +} + +KOKKOS_INLINE_FUNCTION +double AcousticSolver:: +acoustic_dt(const Kokkos::View<const double*>& Vj, + const Kokkos::View<const double*>& cj) const +{ + double dt = std::numeric_limits<double>::max(); + + Kokkos::View<double*> Vj_cj("Vj_cj", m_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); + + return dt; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolver::inverse(const Kokkos::View<const double*>& x) const +{ + Kokkos::View<double*> inv_x("inv_x", x.size()); + Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) { + inv_x(r) = 1./x(r); + }); + + return inv_x; +} + +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); + Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes); +} + +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*> 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",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(m_nr, KOKKOS_LAMBDA(const int& r){ + node_nb_cells(r) = 2; + }); + node_nb_cells(0) = 1; + node_nb_cells(m_nr-1) = 1; + + node_cells(0,0) = 0; + 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(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(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) { + 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 Kokkos::View<const double*> inv_mj=inverse(mj); + + const double tmax=0.2; + double t=0; + + int itermax=std::numeric_limits<int>::max(); + int iteration=0; + + 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; + }); + + 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; + Kokkos::View<double*[2]> Fjr; + + computeExplicitFluxes(xr, xj, + rhoj, uj, pj, cj, Vj, Cjr, + cell_nodes, node_cells, node_nb_cells, node_cell_local_node, + ur, Fjr); + + 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]; + } + uj[j] -= dt*inv_mj[j]*(momentum_fluxes); + Ej[j] -= dt*inv_mj[j]*(energy_fluxes); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { + ej[j] = Ej[j] - 0.5 * uj[j]*uj[j]; + }); + + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ + xr[r] += dt*ur[r]; + }); + + 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::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 new file mode 100644 index 0000000000000000000000000000000000000000..d8429d5acc023acea9ef6930c8de1306f0bc1370 --- /dev/null +++ b/experimental/AcousticSolver.hpp @@ -0,0 +1,75 @@ +#ifndef ACOUSTIC_SOLVER_HPP +#define ACOUSTIC_SOLVER_HPP + +#include <Kokkos_Core.hpp> + +class AcousticSolver +{ +private: + + 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; + + 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 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); + + struct ReduceMin; + + const long int m_nj; + const long int m_nr; + +public: + AcousticSolver(const long int& nj); +}; + +#endif // ACOUSTIC_SOLVER_HPP diff --git a/experimental/AcousticSolverTest.cpp b/experimental/AcousticSolverTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5c1e190b8f06163cb952db133523067d9b142bca --- /dev/null +++ b/experimental/AcousticSolverTest.cpp @@ -0,0 +1,383 @@ +#include <AcousticSolverTest.hpp> +#include <rang.hpp> + +#include <memory> + +#include <BlockPerfectGas.hpp> + +typedef const double my_double; + +struct AcousticSolverTest::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 +const Kokkos::View<const double*> +AcousticSolverTest::computeRhoCj(const Kokkos::View<const double*>& rhoj, + const Kokkos::View<const double*>& cj) +{ + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + m_rhocj[j] = rhoj[j]*cj[j]; + }); + return m_rhocj; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*[2]> +AcousticSolverTest::computeAjr(const Kokkos::View<const double*>& rhocj, + const Kokkos::View<const double*[2]>& Cjr) +{ + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + for (int r=0; r<2; ++r) { + m_Ajr(j,r) = rhocj(j)*Cjr(j,r)*Cjr(j,r); + } + }); + + return m_Ajr; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolverTest::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::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); + const int R = node_cell_local_node(r,j); + sum += Ajr(J,R); + } + m_Ar(r) = sum; + }); + + return m_Ar; +} + +KOKKOS_INLINE_FUNCTION +const Kokkos::View<const double*> +AcousticSolverTest::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::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); + const int R = node_cell_local_node(r,j); + sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj(J); + } + m_br(r) = sum; + }); + + return m_br; +} + +KOKKOS_INLINE_FUNCTION +Kokkos::View<double*> +AcousticSolverTest::computeUr(const Kokkos::View<const double*>& Ar, + const Kokkos::View<const double*>& br) +{ + inverse(Ar, m_inv_Ar); + const Kokkos::View<const double*> invAr = m_inv_Ar; + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) { + m_ur[r]=br(r)*invAr(r); + }); + m_ur[0]=0; + m_ur[m_nr-1]=0; + + return m_ur; +} + +KOKKOS_INLINE_FUNCTION +Kokkos::View<double*[2]> +AcousticSolverTest::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::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { + for (int r=0; r<2; ++r) { + m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+Cjr(j,r)*pj(j); + } + }); + + return m_Fjr; +} + +KOKKOS_INLINE_FUNCTION +double AcousticSolverTest:: +acoustic_dt(const Kokkos::View<const double*>& Vj, + const Kokkos::View<const double*>& cj) const +{ + Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ + m_Vj_over_cj[j] = Vj[j]/cj[j]; + }); + + double dt = std::numeric_limits<double>::max(); + Kokkos::parallel_reduce(m_nj, ReduceMin(m_Vj_over_cj), dt); + + return dt; +} + +KOKKOS_INLINE_FUNCTION +void +AcousticSolverTest::inverse(const Kokkos::View<const double*>& x, + Kokkos::View<double*>& inv_x) const +{ + Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) { + inv_x(r) = 1./x(r); + }); +} + +KOKKOS_INLINE_FUNCTION +void AcousticSolverTest::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); + Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes); +} + +AcousticSolverTest::AcousticSolverTest(const long int& nj) + : m_nj(nj), + m_nr(nj+1), + m_br("br", m_nr), + m_Ajr("Ajr", m_nj), + m_Ar("Ar", m_nr), + m_inv_Ar("inv_Ar", m_nr), + m_Fjr("Fjr", m_nj), + m_ur("ur", m_nr), + m_rhocj("rho_c", m_nj), + m_Vj_over_cj("Vj_over_cj", m_nj) +{ + 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*> 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",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(m_nr, KOKKOS_LAMBDA(const int& r){ + node_nb_cells(r) = 2; + }); + node_nb_cells(0) = 1; + node_nb_cells(m_nr-1) = 1; + + node_cells(0,0) = 0; + 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(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(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) { + 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]; + }); + + Kokkos::View<double*> inv_mj("inv_mj",m_nj); + inverse(mj, inv_mj); + + const double tmax=0.2; + double t=0; + + int itermax=std::numeric_limits<int>::max(); + int iteration=0; + + 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; + }); + + 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; + } + + computeExplicitFluxes(xr, xj, + rhoj, uj, pj, cj, Vj, Cjr, + cell_nodes, node_cells, node_nb_cells, node_cell_local_node, + m_ur, m_Fjr); + + const Kokkos::View<const double*[2]> Fjr = m_Fjr; + const Kokkos::View<const double*> ur = m_ur; + + 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]; + } + uj[j] -= dt*inv_mj[j]*(momentum_fluxes); + Ej[j] -= dt*inv_mj[j]*(energy_fluxes); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) { + ej[j] = Ej[j] - 0.5 * uj[j]*uj[j]; + }); + + Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){ + xr[r] += dt*ur[r]; + }); + + 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::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/AcousticSolverTest.hpp b/experimental/AcousticSolverTest.hpp new file mode 100644 index 0000000000000000000000000000000000000000..f1281635550cb9323452c5d4cb85e688b6216704 --- /dev/null +++ b/experimental/AcousticSolverTest.hpp @@ -0,0 +1,84 @@ +#ifndef ACOUSTIC_SOLVER_TEST_HPP +#define ACOUSTIC_SOLVER_TEST_HPP + +#include <Kokkos_Core.hpp> + +class AcousticSolverTest +{ +private: + 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); + + void inverse(const Kokkos::View<const double*>& x, + Kokkos::View<double*>& inv_x) 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 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); + + struct ReduceMin; + + const long int m_nj; + const long int m_nr; + + Kokkos::View<double*> m_br; + Kokkos::View<double*[2]> m_Ajr; + Kokkos::View<double*> m_Ar; + Kokkos::View<double*> m_inv_Ar; + Kokkos::View<double*[2]> m_Fjr; + Kokkos::View<double*> m_ur; + Kokkos::View<double*> m_rhocj; + Kokkos::View<double*> m_Vj_over_cj; + + +public: + AcousticSolverTest(const long int& nj); +}; + +#endif // ACOUSTIC_SOLVER_TEST_HPP diff --git a/experimental/BlockPerfectGas.hpp b/experimental/BlockPerfectGas.hpp new file mode 100644 index 0000000000000000000000000000000000000000..c018399e5d0644762a2474f3a99da6de20f82ffa --- /dev/null +++ b/experimental/BlockPerfectGas.hpp @@ -0,0 +1,61 @@ +#ifndef BLOCK_PERFECTGAS_HPP +#define BLOCK_PERFECTGAS_HPP + +struct BlockPerfectGas +{ +private: + Kokkos::View<double*> m_rhoj; + Kokkos::View<double*> m_ej; + Kokkos::View<double*> m_pj; + Kokkos::View<double*> m_gammaj; + Kokkos::View<double*> m_cj; + +public: + BlockPerfectGas(Kokkos::View<double*> rhoj, + Kokkos::View<double*> ej, + Kokkos::View<double*> pj, + Kokkos::View<double*> gammaj, + Kokkos::View<double*> cj) + : m_rhoj(rhoj), + m_ej(ej), + m_pj(pj), + m_gammaj(gammaj), + m_cj(cj) + { + ; + } + + void updatePandCFromRhoE() + { + const int nj = m_ej.size(); + const Kokkos::View<const double*> rho = m_rhoj; + const Kokkos::View<const double*> e = m_ej; + const Kokkos::View<const double*> gamma = m_gammaj; + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + const double gamma_minus_one = gamma[j]-1; + m_pj[j] = gamma_minus_one*rho[j]*e[j]; + m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]); + }); + } + + void updateEandCFromRhoP() + { + const int nj = m_ej.size(); + const Kokkos::View<const double*> rho = m_rhoj; + const Kokkos::View<const double*> p = m_pj; + const Kokkos::View<const double*> gamma = m_gammaj; + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + m_ej[j] = p[j]/(rho[j]*(gamma[j]-1)); + }); + + const Kokkos::View<const double*> e = m_ej; + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]); + }); + } +}; + + +#endif // BLOCK_PERFECTGAS_HPP diff --git a/experimental/CMakeLists.txt b/experimental/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..35eeea33d4e73abd2217fa5f8e57f754ec52d3be --- /dev/null +++ b/experimental/CMakeLists.txt @@ -0,0 +1,17 @@ +include_directories(${CMAKE_CURRENT_SOURCE_DIR}) +include_directories(${CMAKE_CURRENT_BINARY_DIR}) + +include_directories(${Kokkos_INCLUDE_DIRS_RET}) +include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include) + +# ------------------- Source files -------------------- + +add_library( + PastisExperimental + AcousticSolver.cpp + AcousticSolverTest.cpp + RawKokkosAcousticSolver.cpp + MeshLessAcousticSolver.cpp + ) + + diff --git a/experimental/MeshLessAcousticSolver.cpp b/experimental/MeshLessAcousticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bf501448130f9e2291a2692d956f2757f5ccfc05 --- /dev/null +++ b/experimental/MeshLessAcousticSolver.cpp @@ -0,0 +1,244 @@ +#include <MeshLessAcousticSolver.hpp> +#include <rang.hpp> + +#include <memory> + +#include <BlockPerfectGas.hpp> + +typedef const double my_double; + +struct MeshLessAcousticSolver::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 MeshLessAcousticSolver:: +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 MeshLessAcousticSolver::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, + Kokkos::View<double*>& ur, + Kokkos::View<double*>& pr) const +{ + // calcul de ur + ur[0]=0; + const size_t nr = ur.size(); + const size_t nj = uj.size(); + + Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) { + const int r = j+1; + const int k = r; + + const double ujr = uj[j]; + const double ukr = uj[k]; + const double pjr = pj[j]; + const double pkr = pj[k]; + + ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]); + }); + 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]); + }); +} + + +MeshLessAcousticSolver::MeshLessAcousticSolver(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::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + xj[j] = 0.5*(xr[j]+xr[j+1]); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + Vj[j] = xr[j+1]-xr[j]; + }); + + 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]; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + inv_mj[j] = 1./mj[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, + 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=j; + int rp=j+1; + + new_uj[j] = uj[j] + dt*inv_mj[j]*(pr[rm]-pr[rp]); + new_Ej[j] = Ej[j] + dt*inv_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/MeshLessAcousticSolver.hpp b/experimental/MeshLessAcousticSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..ccdba58383fa691ed713887a3b0fd09e37b92446 --- /dev/null +++ b/experimental/MeshLessAcousticSolver.hpp @@ -0,0 +1,31 @@ +#ifndef MESH_LESS_ACOUSTIC_SOLVER_HPP +#define MESH_LESS_ACOUSTIC_SOLVER_HPP + +#include <Kokkos_Core.hpp> + +class MeshLessAcousticSolver +{ +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, + Kokkos::View<double*>& ur, + Kokkos::View<double*>& pr) const; + + struct ReduceMin; +public: + MeshLessAcousticSolver(const long int& nj); +}; + +#endif // MESH_LESS_ACOUSTIC_SOLVER_HPP diff --git a/experimental/RawKokkosAcousticSolver.cpp b/experimental/RawKokkosAcousticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d27549d349183d6eb714be477750ca4fee6381b5 --- /dev/null +++ b/experimental/RawKokkosAcousticSolver.cpp @@ -0,0 +1,248 @@ +#include <RawKokkosAcousticSolver.hpp> +#include <Kokkos_Core.hpp> +#include <rang.hpp> + +namespace RawKokkos +{ + +inline double e(double rho, double p, double gamma) +{ + return p/(rho*(gamma-1)); +} + +inline double p(double rho, double e, double gamma) +{ + return (gamma-1)*rho*e; +} + +typedef const double my_double; + +struct 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(); + } +}; + + +double acoustic_dt(const Kokkos::View<double*>& Vj, + const Kokkos::View<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::parallel_reduce(nj, ReduceMin(Vj_cj), dt); + + // Kokkos::parallel_reduce(n, KOKKOS_LAMBDA(const long i, long& lcount) { + // lcount += (i % 2) == 0; + // }, count); + return dt; +} + + +void computeExplicitFluxes(const Kokkos::View<double*>& xr, + const Kokkos::View<double*>& xj, + const Kokkos::View<double*>& rhoj, + const Kokkos::View<double*>& uj, + const Kokkos::View<double*>& pj, + const Kokkos::View<double*>& cj, + const Kokkos::View<double*>& Vj, + Kokkos::View<double*>& ur, + Kokkos::View<double*>& pr) +{ + // calcul de ur + ur[0]=0; + const size_t nr = ur.size(); + const size_t nj = uj.size(); + + Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) { + const int r = j+1; + const int k = r; + const double ujr = uj[j]; + const double ukr = uj[k]; + const double pjr = pj[j]; + const double pkr = pj[k]; + + ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]); + }); + 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]); + }); +} + + +void 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::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + xj[j] = 0.5*(xr[j]+xr[j+1]); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + Vj[j] = xr[j+1]-xr[j]; + }); + + 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; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + ej[j] = e(rhoj[j],pj[j],gammaj[j]); + }); + + 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){ + cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + mj[j] = rhoj[j] * Vj[j]; + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + inv_mj[j] = 1./mj[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, + ur, pr); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + int rm=j; + int rp=j+1; + + uj[j] += dt*inv_mj[j]*(pr[rm]-pr[rp]); + Ej[j] += dt*inv_mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]); + }); + + Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ + ej[j] = Ej[j] - 0.5 * uj[j]*uj[j]; + }); + + Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ + xr[r] += dt*ur[r]; + }); + + 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]; + pj[j] = p(rhoj[j], ej[j], gammaj[j]); + cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); // inv_mj*vj + }); + + ++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/RawKokkosAcousticSolver.hpp b/experimental/RawKokkosAcousticSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e9c36ff437a7ef0f0c28c72e93b0a367d8142185 --- /dev/null +++ b/experimental/RawKokkosAcousticSolver.hpp @@ -0,0 +1,9 @@ +#ifndef RAW_KOKKOS_ACOUSTIC_SOLVER_HPP +#define RAW_KOKKOS_ACOUSTIC_SOLVER_HPP + +namespace RawKokkos +{ + void AcousticSolver(const long int& nj); +} + +#endif // RAW_KOKKOS_ACOUSTIC_SOLVER_HPP diff --git a/main.cpp b/main.cpp index fc1e3316a0b6073b51232fffa897e406d706380c..fb203f33487c79a86fdaef0480cdfc3dd7b42585 100644 --- a/main.cpp +++ b/main.cpp @@ -6,120 +6,15 @@ #include <SignalManager.hpp> #include <ConsoleManager.hpp> +#include <RawKokkosAcousticSolver.hpp> +#include <MeshLessAcousticSolver.hpp> +#include <AcousticSolver.hpp> +#include <AcousticSolverTest.hpp> + #include <CLI/CLI.hpp> #include <cassert> #include <limits> - -inline double e(double rho, double p, double gamma) -{ - return p/(rho*(gamma-1)); -} - -inline double p(double rho, double e, double gamma) -{ - return (gamma-1)*rho*e; -} - -typedef const double my_double; - -struct 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(); - } -}; - - -double acoustic_dt(const Kokkos::View<double*>& Vj, - const Kokkos::View<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::parallel_reduce(nj, ReduceMin(Vj_cj), dt); - - // Kokkos::parallel_reduce(n, KOKKOS_LAMBDA(const long i, long& lcount) { - // lcount += (i % 2) == 0; - // }, count); - return dt; -} - - -void computeExplicitFluxes(const Kokkos::View<double*>& xr, - const Kokkos::View<double*>& xj, - const Kokkos::View<double*>& rhoj, - const Kokkos::View<double*>& uj, - const Kokkos::View<double*>& pj, - const Kokkos::View<double*>& cj, - const Kokkos::View<double*>& Vj, - Kokkos::View<double*>& ur, - Kokkos::View<double*>& pr) -{ - // calcul de ur - ur[0]=0; - const size_t nr = ur.size(); - const size_t nj = uj.size(); - - Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) { - const int r = j+1; - const int k = r; - const double ujr = uj[j]; - const double ukr = uj[k]; - const double pjr = pj[j]; - const double pkr = pj[k]; - - ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]); - }); - 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]); - }); -} - - -#warning clean-up and add warning message when release version is run +#include <map> int main(int argc, char *argv[]) { @@ -182,148 +77,53 @@ int main(int argc, char *argv[]) Kokkos::initialize(argc, argv); Kokkos::DefaultExecutionSpace::print_configuration(std::cout); - const long& nj=number; - - 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::Timer timer; - timer.reset(); - - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ - xr[r] = r*delta_x; - }); - - - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - xj[j] = 0.5*(xr[j]+xr[j+1]); - }); - - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - Vj[j] = xr[j+1]-xr[j]; - }); - - 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; - }); - - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - ej[j] = e(rhoj[j],pj[j],gammaj[j]); - }); + std::map<std::string, double> method_cost_map; - 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){ - cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[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; - } - - if (iteration%100 == 0) { - std::cout << "dt=" << dt << " t=" << t << " i=" << iteration << std::endl; - } - - Kokkos::View<double*> ur("ur", nr); - Kokkos::View<double*> pr("pr", nr); + { // Basic function based acoustic solver + Kokkos::Timer timer; + timer.reset(); + RawKokkos::AcousticSolver(number); + method_cost_map["RawKokkos"] = timer.seconds(); + } - computeExplicitFluxes(xr, xj, - rhoj, uj, pj, cj, Vj, - ur, pr); - - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - int rm=j; - int rp=j+1; + { // class for acoustic solver (mesh less) + Kokkos::Timer timer; + timer.reset(); + MeshLessAcousticSolver acoustic_solver(number); + method_cost_map["MeshLessAcousticSolver"] = timer.seconds(); + } - uj[j] += dt/mj[j]*(pr[rm]-pr[rp]); - Ej[j] += dt/mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]); - }); + { // class for acoustic solver test + Kokkos::Timer timer; + timer.reset(); + AcousticSolverTest acoustic_solver(number); + method_cost_map["AcousticSolverTest"] = timer.seconds(); + } - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - ej[j] = Ej[j] - 0.5 * uj[j]*uj[j]; - }); + { // class for acoustic solver + Kokkos::Timer timer; + timer.reset(); + AcousticSolver acoustic_solver(number); + method_cost_map["AcousticSolver"] = timer.seconds(); + } - Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){ - xr[r] += dt*ur[r]; - }); + Kokkos::finalize(); - 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]; - }); + std::cout << "----------------------\n"; - Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){ - rhoj[j] = mj[j]/Vj[j]; - pj[j] = p(rhoj[j], ej[j], gammaj[j]); - cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); // inv_mj*vj - }); - - ++iteration; + std::string::size_type size=0; + for (const auto& method_cost : method_cost_map) { + size = std::max(size, method_cost.first.size()); } - std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset - << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; - double count_time = timer.seconds(); - std::cout << "* Execution time: " << rang::style::bold << count_time << rang::style::reset << '\n'; - - { - std::ofstream fout("rho"); - - for (int j=0; j<nj; ++j) { - fout << xj[j] << ' ' << rhoj[j] << '\n'; - } + for (const auto& method_cost : method_cost_map) { + std::cout << "* [" << std::setw(size) << std::left + << method_cost.first + << " ] Execution time: " + << rang::style::bold + << method_cost.second + << rang::style::reset << '\n'; } - Kokkos::finalize(); return 0; }