Skip to content
Snippets Groups Projects
Commit 9ee1d44d authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Began acoustic solver with explicit connectivity

parent 7a8e79f6
No related branches found
No related tags found
No related merge requests found
#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";
}
#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
......@@ -8,6 +8,7 @@ include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include)
add_library(
PastisExperimental
AcousticSolver.cpp
RawKokkosAcousticSolver.cpp
MeshLessAcousticSolver.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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment