Skip to content

Commit

Permalink
Merge pull request #7 from pkestene/fix/mhd-cfl-bug
Browse files Browse the repository at this point in the history 
Fix cfl computation for MHD.
  • Loading branch information
@pkestene
pkestene authored Nov 8, 2024
2 parents da92099 + e702132 commit 5f0c8c7
Show file tree
Hide file tree
Showing 4 changed files with 50 additions and 52 deletions.
4 changes: 2 additions & 2 deletions src/muscl/MHDRunFunctors2D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,9 @@ class ComputeDtFunctor2D_MHD : public MHDBaseFunctor2D

// static method which does it all: create and execute functor
static void
apply(HydroParams params, DataArray2d Udata, real_t & invDt)
apply(HydroParams params, DataArray2d Qdata, real_t & invDt)
{
ComputeDtFunctor2D_MHD functor(params, Udata);
ComputeDtFunctor2D_MHD functor(params, Qdata);
Kokkos::Max<real_t> reducer(invDt);
Kokkos::parallel_reduce(
"ComputeDtFunctor2D_MHD",
Expand Down
4 changes: 2 additions & 2 deletions src/muscl/MHDRunFunctors3D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,9 @@ class ComputeDtFunctor3D_MHD : public MHDBaseFunctor3D

// static method which does it all: create and execute functor
static void
apply(HydroParams params, DataArray3d Udata, real_t & invDt)
apply(HydroParams params, DataArray3d Qdata, real_t & invDt)
{
ComputeDtFunctor3D_MHD functor(params, Udata);
ComputeDtFunctor3D_MHD functor(params, Qdata);
Kokkos::Max<real_t> reducer(invDt);
Kokkos::parallel_reduce(Kokkos::MDRangePolicy<Kokkos::Rank<3>>(
{ 0, 0, 0 }, { params.isize, params.jsize, params.ksize }),
Expand Down
55 changes: 31 additions & 24 deletions src/muscl/SolverMHDMuscl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -368,15 +368,9 @@ SolverMHDMuscl<3>::computeEmfAndUpdate(real_t dt, DataArray Udata)
// ///////////////////////////////////////////
template <>
void
SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt)
SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out)
{

real_t dtdx;
real_t dtdy;

dtdx = dt / params.dx;
dtdy = dt / params.dy;

// fill ghost cell in data_in
timers[TIMER_BOUNDARIES]->start();
make_boundaries(data_in);
Expand All @@ -392,17 +386,25 @@ SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
// convert conservative variable into primitives ones for the entire domain
convertToPrimitives(data_in);

// compute new dt
timers[TIMER_DT]->start();
compute_dt();
timers[TIMER_DT]->stop();

const auto dtdx = m_dt / params.dx;
const auto dtdy = m_dt / params.dy;

if (params.implementationVersion == 0)
{

// trace computation: fill arrays qm_x, qm_y, qp_x, qp_y
computeTrace(data_in, dt);
computeTrace(data_in, m_dt);

// Compute flux via Riemann solver and update (time integration)
computeFluxesAndStore(dt);
computeFluxesAndStore(m_dt);

// Compute Emf
computeEmfAndStore(dt);
computeEmfAndStore(m_dt);

// actual update with fluxes
UpdateFunctor2D_MHD::apply(params, data_out, Fluxes_x, Fluxes_y, dtdx, dtdy);
Expand All @@ -413,13 +415,13 @@ SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
else if (params.implementationVersion == 1)
{
// trace computation: fill arrays qm_x, qm_y, qp_x, qp_y
computeTrace(data_in, dt);
computeTrace(data_in, m_dt);

// Compute flux via Riemann solver and update (time integration)
computeFluxesAndUpdate(dt, data_out);
computeFluxesAndUpdate(m_dt, data_out);

// Compute Emf and update magnetic field
computeEmfAndUpdate(dt, data_out);
computeEmfAndUpdate(m_dt, data_out);
}
else if (params.implementationVersion == 2)
{
Expand Down Expand Up @@ -458,13 +460,9 @@ SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
// ///////////////////////////////////////////
template <>
void
SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt)
SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out)
{

const real_t dtdx = dt / params.dx;
const real_t dtdy = dt / params.dy;
const real_t dtdz = dt / params.dz;

// fill ghost cell in data_in
timers[TIMER_BOUNDARIES]->start();
make_boundaries(data_in);
Expand All @@ -480,6 +478,15 @@ SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
// convert conservative variable into primitives ones for the entire domain
convertToPrimitives(data_in);

// compute new dt
timers[TIMER_DT]->start();
compute_dt();
timers[TIMER_DT]->stop();

const auto dtdx = m_dt / params.dx;
const auto dtdy = m_dt / params.dy;
const auto dtdz = m_dt / params.dz;

if (params.implementationVersion == 0)
{

Expand All @@ -490,13 +497,13 @@ SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
computeMagSlopes(data_in);

// trace computation: fill arrays qm_x, qm_y, qm_z, qp_x, qp_y, qp_z
computeTrace(data_in, dt);
computeTrace(data_in, m_dt);

// Compute flux via Riemann solver and update (time integration)
computeFluxesAndStore(dt);
computeFluxesAndStore(m_dt);

// Compute Emf
computeEmfAndStore(dt);
computeEmfAndStore(m_dt);

// actual update with fluxes
UpdateFunctor3D_MHD::apply(params, data_out, Fluxes_x, Fluxes_y, Fluxes_z, dtdx, dtdy, dtdz);
Expand All @@ -513,13 +520,13 @@ SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, r
computeMagSlopes(data_in);

// trace computation: fill arrays qm_x, qm_y, qm_z, qp_x, qp_y, qp_z
computeTrace(data_in, dt);
computeTrace(data_in, m_dt);

// Compute flux via Riemann solver and update (time integration)
computeFluxesAndUpdate(dt, data_out);
computeFluxesAndUpdate(m_dt, data_out);

// Compute Emf and update magnetic field
computeEmfAndUpdate(dt, data_out);
computeEmfAndUpdate(m_dt, data_out);
}
else if (params.implementationVersion == 2)
{
Expand Down
39 changes: 15 additions & 24 deletions src/muscl/SolverMHDMuscl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,10 +157,10 @@ class SolverMHDMuscl : public euler_kokkos::SolverBase

//! numerical scheme
void
godunov_unsplit(real_t dt);
godunov_unsplit();

void
godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt);
godunov_unsplit_impl(DataArray data_in, DataArray data_out);

void
convertToPrimitives(DataArray Udata);
Expand Down Expand Up @@ -390,6 +390,9 @@ SolverMHDMuscl<dim>::SolverMHDMuscl(HydroParams & params, ConfigMap & configMap)
// copy U into U2
Kokkos::deep_copy(U2, U);

// primitive variables are necessary for computing time step
convertToPrimitives(U);

// compute initialize time step
compute_dt();

Expand Down Expand Up @@ -717,22 +720,15 @@ double
SolverMHDMuscl<dim>::compute_dt_local()
{

real_t dt;
real_t invDt = ZERO_F;
DataArray Udata;

// which array is the current one ?
if (m_iteration % 2 == 0)
Udata = U;
else
Udata = U2;
real_t dt;
real_t invDt = ZERO_F;

// alias to actual device functor
using ComputeDtFunctor =
typename std::conditional<dim == 2, ComputeDtFunctor2D_MHD, ComputeDtFunctor3D_MHD>::type;

// call device functor
ComputeDtFunctor::apply(params, Udata, invDt);
ComputeDtFunctor::apply(params, Q, invDt);

dt = params.settings.cfl / invDt;

Expand Down Expand Up @@ -778,13 +774,8 @@ SolverMHDMuscl<dim>::next_iteration_impl()
} // end output
} // end enable output

// compute new dt
timers[TIMER_DT]->start();
compute_dt();
timers[TIMER_DT]->stop();

// perform one step integration
godunov_unsplit(m_dt);
godunov_unsplit();

} // SolverMHDMuscl::next_iteration_impl

Expand All @@ -795,16 +786,16 @@ SolverMHDMuscl<dim>::next_iteration_impl()
// ///////////////////////////////////////////
template <int dim>
void
SolverMHDMuscl<dim>::godunov_unsplit(real_t dt)
SolverMHDMuscl<dim>::godunov_unsplit()
{

if (m_iteration % 2 == 0)
{
godunov_unsplit_impl(U, U2, dt);
godunov_unsplit_impl(U, U2);
}
else
{
godunov_unsplit_impl(U2, U, dt);
godunov_unsplit_impl(U2, U);
}

} // SolverMHDMuscl::godunov_unsplit
Expand All @@ -816,7 +807,7 @@ SolverMHDMuscl<dim>::godunov_unsplit(real_t dt)
// ///////////////////////////////////////////
template <int dim>
void
SolverMHDMuscl<dim>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt)
SolverMHDMuscl<dim>::godunov_unsplit_impl(DataArray data_in, DataArray data_out)
{

// 2d / 3d implementation are specialized
Expand All @@ -826,12 +817,12 @@ SolverMHDMuscl<dim>::godunov_unsplit_impl(DataArray data_in, DataArray data_out,
// 2d
template <>
void
SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt);
SolverMHDMuscl<2>::godunov_unsplit_impl(DataArray data_in, DataArray data_out);

// 3d
template <>
void
SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt);
SolverMHDMuscl<3>::godunov_unsplit_impl(DataArray data_in, DataArray data_out);


// =======================================================
Expand Down

0 comments on commit 5f0c8c7

Please sign in to comment.