Merge pull request #3 from pkestene/fix/mhd-cfl-bug

pkestene · web-flow · commit c5fa7dcb5f9b · 2024-11-09T14:09:45.000+01:00
Fix cfl computation for MHD.
diff --git 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<real_t>    reducer(invDt);
     Kokkos::parallel_reduce(
       "ComputeDtFunctor2D_MHD",
diff --git 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<real_t>    reducer(invDt);
     Kokkos::parallel_reduce(Kokkos::MDRangePolicy<Kokkos::Rank<3>>(
                               { 0, 0, 0 }, { params.isize, params.jsize, params.ksize }),
diff --git a/src/muscl/SolverHydroMuscl.cpp b/src/muscl/SolverHydroMuscl.cpp
@@ -2,9 +2,6 @@
 #include <cstdio>
 #include <cstdbool>
 #include <iostream> // for std::cout
-#include <sstream>
-#include <fstream>
-#include <algorithm>
 
 #include "muscl/SolverHydroMuscl.h"
 #include "shared/HydroParams.h"
diff --git a/src/muscl/SolverMHDMuscl.cpp b/src/muscl/SolverMHDMuscl.cpp
@@ -372,15 +372,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);
@@ -396,17 +390,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);
@@ -417,13 +419,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)
   {
@@ -462,13 +464,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);
@@ -484,6 +482,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)
   {
 
@@ -494,13 +501,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);
@@ -517,13 +524,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
@@ -10,9 +10,6 @@
 #include <cstdio>
 #include <cstdbool>
 #include <iostream>
-#include <sstream>
-#include <fstream>
-#include <algorithm>
 
 // shared
 #include "shared/SolverBase.h"
@@ -165,10 +162,10 @@ class SolverMHDMuscl : public ppkMHD::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);
@@ -398,6 +395,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();
 
@@ -724,22 +724,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;
 
@@ -785,13 +778,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
 
@@ -802,16 +790,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
@@ -823,7 +811,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
@@ -833,12 +821,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);
 
 
 // =======================================================
