diff --git a/src/muscl/MHDRunFunctors2D.h b/src/muscl/MHDRunFunctors2D.h index 442c13f..5c10702 100644 --- a/src/muscl/MHDRunFunctors2D.h +++ b/src/muscl/MHDRunFunctors2D.h @@ -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 reducer(invDt); Kokkos::parallel_reduce( "ComputeDtFunctor2D_MHD", diff --git a/src/muscl/MHDRunFunctors3D.h b/src/muscl/MHDRunFunctors3D.h index d119923..5c995cb 100644 --- a/src/muscl/MHDRunFunctors3D.h +++ b/src/muscl/MHDRunFunctors3D.h @@ -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 reducer(invDt); Kokkos::parallel_reduce(Kokkos::MDRangePolicy>( { 0, 0, 0 }, { params.isize, params.jsize, params.ksize }), diff --git a/src/muscl/SolverMHDMuscl.cpp b/src/muscl/SolverMHDMuscl.cpp index fe2d42c..4d98859 100644 --- a/src/muscl/SolverMHDMuscl.cpp +++ b/src/muscl/SolverMHDMuscl.cpp @@ -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); @@ -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); @@ -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) { @@ -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); @@ -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) { @@ -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); @@ -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) { diff --git a/src/muscl/SolverMHDMuscl.h b/src/muscl/SolverMHDMuscl.h index 08966e4..4d5f1c0 100644 --- a/src/muscl/SolverMHDMuscl.h +++ b/src/muscl/SolverMHDMuscl.h @@ -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); @@ -390,6 +390,9 @@ SolverMHDMuscl::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(); @@ -717,22 +720,15 @@ double SolverMHDMuscl::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::type; // call device functor - ComputeDtFunctor::apply(params, Udata, invDt); + ComputeDtFunctor::apply(params, Q, invDt); dt = params.settings.cfl / invDt; @@ -778,13 +774,8 @@ SolverMHDMuscl::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 @@ -795,16 +786,16 @@ SolverMHDMuscl::next_iteration_impl() // /////////////////////////////////////////// template void -SolverMHDMuscl::godunov_unsplit(real_t dt) +SolverMHDMuscl::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 @@ -816,7 +807,7 @@ SolverMHDMuscl::godunov_unsplit(real_t dt) // /////////////////////////////////////////// template void -SolverMHDMuscl::godunov_unsplit_impl(DataArray data_in, DataArray data_out, real_t dt) +SolverMHDMuscl::godunov_unsplit_impl(DataArray data_in, DataArray data_out) { // 2d / 3d implementation are specialized @@ -826,12 +817,12 @@ SolverMHDMuscl::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); // =======================================================