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

externalizes Euler integration from acoustic solver

parent 1decef79
Branches
Tags
No related merge requests found
...@@ -182,19 +182,6 @@ private: ...@@ -182,19 +182,6 @@ private:
}); });
} }
KOKKOS_INLINE_FUNCTION
double acoustic_dt(const Kokkos::View<const double*>& Vj,
const Kokkos::View<const double*>& cj) const {
Kokkos::parallel_for(m_mesh.numberOfCells(), 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_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt);
return dt;
}
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr, void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
const Kokkos::View<const Rd*>& xj, const Kokkos::View<const Rd*>& xj,
...@@ -225,7 +212,6 @@ private: ...@@ -225,7 +212,6 @@ private:
Kokkos::View<double*> m_rhocj; Kokkos::View<double*> m_rhocj;
Kokkos::View<double*> m_Vj_over_cj; Kokkos::View<double*> m_Vj_over_cj;
public: public:
AcousticSolverWithMesh(MeshData& mesh_data, AcousticSolverWithMesh(MeshData& mesh_data,
UnknownsType& unknowns) UnknownsType& unknowns)
...@@ -240,6 +226,26 @@ public: ...@@ -240,6 +226,26 @@ public:
m_ur("ur", m_mesh.numberOfNodes()), m_ur("ur", m_mesh.numberOfNodes()),
m_rhocj("rho_c", m_mesh.numberOfCells()), m_rhocj("rho_c", m_mesh.numberOfCells()),
m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells()) m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
{
;
}
KOKKOS_INLINE_FUNCTION
double acoustic_dt(const Kokkos::View<const double*>& Vj,
const Kokkos::View<const double*>& cj) const {
Kokkos::parallel_for(m_mesh.numberOfCells(), 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_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt);
return dt;
}
void computeNextStep(const double& t, const double& dt,
UnknownsType& unknowns)
{ {
Kokkos::View<double*> rhoj = unknowns.rhoj(); Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<Rd*> uj = unknowns.uj(); Kokkos::View<Rd*> uj = unknowns.uj();
...@@ -249,35 +255,11 @@ public: ...@@ -249,35 +255,11 @@ public:
Kokkos::View<double*> pj = unknowns.pj(); Kokkos::View<double*> pj = unknowns.pj();
Kokkos::View<double*> gammaj = unknowns.gammaj(); Kokkos::View<double*> gammaj = unknowns.gammaj();
Kokkos::View<double*> cj = unknowns.cj(); Kokkos::View<double*> cj = unknowns.cj();
const Kokkos::View<const double*> mj = unknowns.mj();
Kokkos::View<Rd*> xr = m_mesh.xr();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
Kokkos::View<Rd*> xr = m_mesh.xr();
Kokkos::View<double*> inv_mj("inv_mj",m_mesh.numberOfCells());
inverse(mj, inv_mj);
const double tmax=0.2;
double t=0;
int itermax=std::numeric_limits<int>::max();
int iteration=0;
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
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); computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr);
...@@ -285,7 +267,10 @@ public: ...@@ -285,7 +267,10 @@ public:
const Kokkos::View<const Rd*> ur = m_ur; const Kokkos::View<const Rd*> ur = m_ur;
const Kokkos::View<const unsigned short*> cell_nb_nodes const Kokkos::View<const unsigned short*> cell_nb_nodes
= m_connectivity.cellNbNodes(); = m_connectivity.cellNbNodes();
const Kokkos::View<const unsigned int**>& cell_nodes
= m_connectivity.cellNodes();
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
Rd momentum_fluxes = zero; Rd momentum_fluxes = zero;
double energy_fluxes = 0; double energy_fluxes = 0;
...@@ -308,18 +293,10 @@ public: ...@@ -308,18 +293,10 @@ public:
m_mesh_data.updateAllData(); m_mesh_data.updateAllData();
const Kokkos::View<const double*> mj = unknowns.mj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
rhoj[j] = mj[j]/Vj[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";
} }
}; };
......
...@@ -24,6 +24,7 @@ private: ...@@ -24,6 +24,7 @@ private:
Kokkos::View<double*> m_gammaj; Kokkos::View<double*> m_gammaj;
Kokkos::View<double*> m_cj; Kokkos::View<double*> m_cj;
Kokkos::View<double*> m_mj; Kokkos::View<double*> m_mj;
Kokkos::View<double*> m_inv_mj;
public: public:
Kokkos::View<double*> rhoj() Kokkos::View<double*> rhoj()
...@@ -106,6 +107,16 @@ public: ...@@ -106,6 +107,16 @@ public:
return m_mj; return m_mj;
} }
Kokkos::View<double*> invMj()
{
return m_inv_mj;
}
const Kokkos::View<const double*> invMj() const
{
return m_inv_mj;
}
void initializeSod() void initializeSod()
{ {
const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
...@@ -146,6 +157,9 @@ public: ...@@ -146,6 +157,9 @@ public:
m_mj[j] = m_rhoj[j] * Vj[j]; m_mj[j] = m_rhoj[j] * Vj[j];
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_inv_mj[j] = 1./m_mj[j];
});
} }
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data) FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
...@@ -158,7 +172,8 @@ public: ...@@ -158,7 +172,8 @@ public:
m_pj("pj",m_mesh.numberOfCells()), m_pj("pj",m_mesh.numberOfCells()),
m_gammaj("gammaj",m_mesh.numberOfCells()), m_gammaj("gammaj",m_mesh.numberOfCells()),
m_cj("cj",m_mesh.numberOfCells()), m_cj("cj",m_mesh.numberOfCells()),
m_mj("mj",m_mesh.numberOfCells()) m_mj("mj",m_mesh.numberOfCells()),
m_inv_mj("inv_mj",m_mesh.numberOfCells())
{ {
; ;
} }
......
...@@ -130,6 +130,41 @@ int main(int argc, char *argv[]) ...@@ -130,6 +130,41 @@ int main(int argc, char *argv[])
AcousticSolverWithMesh<MeshDataType> acoustic_solver(mesh_data, unknowns); AcousticSolverWithMesh<MeshDataType> acoustic_solver(mesh_data, unknowns);
const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.2;
double t=0;
int itermax=std::numeric_limits<int>::max();
int iteration=0;
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<double*> ej = unknowns.ej();
Kokkos::View<double*> pj = unknowns.pj();
Kokkos::View<double*> gammaj = unknowns.gammaj();
Kokkos::View<double*> cj = unknowns.cj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
while((t<tmax) and (iteration<itermax)) {
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
t=tmax;
}
acoustic_solver.computeNextStep(t,dt, unknowns);
block_eos.updatePandCFromRhoE();
t += dt;
++iteration;
}
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
{ // gnuplot output for density { // gnuplot output for density
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment