Add prolongate prim for hydro

Author: pgrete

src/eos/adiabatic_glmmhd.hpp

@@ -15,6 +15,7 @@
 // Athena headers
 #include "../main.hpp"
 #include "eos.hpp"
+#include "utils/error_checking.hpp"
 
 using parthenon::MeshBlock;
 using parthenon::MeshBlockData;
@@ -30,6 +31,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
         gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
+  void PrimitiveToConserved(MeshData<Real> *md) const override {
+    PARTHENON_FAIL("needs impl.");
+  }
 
   KOKKOS_INLINE_FUNCTION
   Real GetGamma() const { return gamma_; }

src/eos/adiabatic_hydro.cpp

@@ -69,3 +69,35 @@ void AdiabaticHydroEOS::ConservedToPrimitive(MeshData<Real> *md) const {
   pkg->UpdateParam<std::int64_t>("fixed_num_cells_floor_temp",
                                  floor_temp_pkg + floor_temp);
 }
+
+//----------------------------------------------------------------------------------------
+// \!fn void EquationOfState::PrimitiveToConserved(
+//           Container<Real> &rc,
+//           int il, int iu, int jl, int ju, int kl, int ku)
+// \brief Converts primitive to conserved variables in adiabatic hydro.
+// Not using any floors here and failing loudly because those fixes
+// are applied in ConsToPrim call. Should discss advantags and disadvantages of this
+// approach.
+void AdiabaticHydroEOS::PrimitiveToConserved(MeshData<Real> *md) const {
+  auto cons_pack = md->PackVariables(std::vector<std::string>{"cons"});
+  auto const prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+  auto ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::entire);
+  auto jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::entire);
+  auto kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::entire);
+
+  auto pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+  const auto nhydro = pkg->Param<int>("nhydro");
+  const auto nscalars = pkg->Param<int>("nscalars");
+
+  auto this_on_device = (*this);
+
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "PrimitiveToConserved", parthenon::DevExecSpace(), 0,
+      cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        auto &cons = cons_pack(b);
+        const auto &prim = prim_pack(b);
+
+        this_on_device.PrimToCons(cons, prim, nhydro, nscalars, k, j, i);
+      });
+}

src/eos/adiabatic_hydro.hpp

@@ -31,6 +31,7 @@ class AdiabaticHydroEOS : public EquationOfState {
         gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
+  void PrimitiveToConserved(MeshData<Real> *md) const override;
 
   KOKKOS_INLINE_FUNCTION
   Real GetGamma() const { return gamma_; }
@@ -151,6 +152,55 @@ class AdiabaticHydroEOS : public EquationOfState {
     return floors_used;
   }
 
+  //----------------------------------------------------------------------------------------
+  // \!fn Real EquationOfState::PrimToCons(View4D cons, View4D prim, const int& k, const
+  // int& j, const int& i) \brief Fills an array of conservatives given an array of
+  // primities, currently without floors
+  template <typename View4D>
+  KOKKOS_INLINE_FUNCTION void PrimToCons(View4D cons, View4D prim, const int &nhydro,
+                                         const int &nscalars, const int &k, const int &j,
+                                         const int &i) const {
+    Real gm1 = GetGamma() - 1.0;
+
+    Real &u_d = cons(IDN, k, j, i);
+    Real &u_m1 = cons(IM1, k, j, i);
+    Real &u_m2 = cons(IM2, k, j, i);
+    Real &u_m3 = cons(IM3, k, j, i);
+    Real &u_e = cons(IEN, k, j, i);
+
+    Real &w_d = prim(IDN, k, j, i);
+    Real &w_vx = prim(IV1, k, j, i);
+    Real &w_vy = prim(IV2, k, j, i);
+    Real &w_vz = prim(IV3, k, j, i);
+    Real &w_p = prim(IPR, k, j, i);
+
+    PARTHENON_REQUIRE(w_d != 0.0,
+                      "Densities should never be exactly 0! This points to working with "
+                      "some default initialized and/or uninitialized data.");
+    PARTHENON_REQUIRE(w_d > 0.0,
+                      "Got negative density. Might need to impl floors here too.");
+    PARTHENON_REQUIRE(w_p != 0.0,
+                      "Pressure should never be exactly 0! This points to working with "
+                      "some default initialized and/or uninitialized data.");
+    PARTHENON_REQUIRE(w_p > 0.0,
+                      "Got negative pressure. Might need to impl floors here too.");
+    // apply density floor, without changing momentum or energy
+    u_d = w_d;
+
+    Real di = 1.0 / u_d;
+    u_m1 = w_d * w_vx;
+    u_m2 = w_d * w_vy;
+    u_m3 = w_d * w_vz;
+
+    const Real e_k = 0.5 * w_d * (SQR(w_vx) + SQR(w_vy) + SQR(w_vz));
+    u_e = e_k + w_p / gm1;
+
+    // Convert passive scalars
+    for (auto n = nhydro; n < nhydro + nscalars; ++n) {
+      cons(n, k, j, i) = prim(n, k, j, i) * w_d;
+    }
+  }
+
  private:
   Real gamma_; // ratio of specific heats
 };

src/eos/eos.hpp

@@ -38,6 +38,7 @@ class EquationOfState {
         internal_e_floor_(internal_e_floor), velocity_ceiling_(velocity_ceiling),
         internal_e_ceiling_(internal_e_ceiling) {}
   virtual void ConservedToPrimitive(MeshData<Real> *md) const = 0;
+  virtual void PrimitiveToConserved(MeshData<Real> *md) const = 0;
 
   KOKKOS_INLINE_FUNCTION
   Real GetPressureFloor() const { return pressure_floor_; }

src/hydro/hydro.cpp

@@ -32,6 +32,7 @@
 #include "diffusion/diffusion.hpp"
 #include "glmmhd/glmmhd.hpp"
 #include "hydro.hpp"
+#include "interface/metadata.hpp"
 #include "interface/params.hpp"
 #include "outputs/outputs.hpp"
 #include "prolongation/custom_ops.hpp"
@@ -217,6 +218,12 @@ void ConsToPrim(MeshData<Real> *md) {
       md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro")->Param<T>("eos");
   eos.ConservedToPrimitive(md);
 }
+template <class T>
+void PrimToCons(MeshData<Real> *md) {
+  const auto &eos =
+      md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro")->Param<T>("eos");
+  eos.PrimitiveToConserved(md);
+}
 
 // Add unsplit sources, i.e., source that are integrated in all stages of the
 // explicit integration scheme.
@@ -493,6 +500,9 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     Units units(pin, pkg);
   }
 
+  const auto prolong_prims = pin->GetOrAddBoolean("hydro", "prolongate_prims", false);
+  pkg->AddParam<>("prolongate_prims", prolong_prims);
+
   auto eos_str = pin->GetString("hydro", "eos");
   if (eos_str == "adiabatic") {
     Real gamma = pin->GetReal("hydro", "gamma");
@@ -717,7 +727,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     if (fluid == Fluid::euler) {
       AdiabaticHydroEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
       pkg->AddParam<>("eos", eos);
-      pkg->FillDerivedMesh = ConsToPrim<AdiabaticHydroEOS>;
+      if (prolong_prims) {
+        pkg->PreCommFillDerivedMesh = ConsToPrim<AdiabaticHydroEOS>;
+        pkg->FillDerivedMesh = PrimToCons<AdiabaticHydroEOS>;
+      } else {
+        pkg->FillDerivedMesh = ConsToPrim<AdiabaticHydroEOS>;
+      }
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::euler>;
     } else if (fluid == Fluid::glmmhd) {
       AdiabaticGLMMHDEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
@@ -787,15 +802,32 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     prim_labels.emplace_back("scalar_" + std::to_string(i));
   }
 
-  Metadata m(
-      {Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes},
-      std::vector<int>({nhydro + nscalars}), cons_labels);
+  Metadata m;
+  // In principle, it'd be nicer to work with .Set(), but see
+  // https://github.com/parthenon-hpc-lab/parthenon/issues/844 so let's
+  // be safe than sorry.
+  if (prolong_prims) {
+    // no FillGhost here
+    m = Metadata({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes},
+                 std::vector<int>({nhydro + nscalars}), cons_labels);
+  } else {
+    m = Metadata({Metadata::Cell, Metadata::Independent, Metadata::FillGhost,
+                  Metadata::WithFluxes},
+                 std::vector<int>({nhydro + nscalars}), cons_labels);
+  }
   m.RegisterRefinementOps<refinement_ops::ProlongateCellMinModMultiD,
                           parthenon::refinement_ops::RestrictAverage>();
   pkg->AddField("cons", m);
 
-  m = Metadata({Metadata::Cell, Metadata::Derived}, std::vector<int>({nhydro + nscalars}),
-               prim_labels);
+  if (prolong_prims) {
+    m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::FillGhost},
+                 std::vector<int>({nhydro + nscalars}), prim_labels);
+    m.RegisterRefinementOps<refinement_ops::ProlongateCellMinModMultiD,
+                            parthenon::refinement_ops::RestrictAverage>();
+  } else {
+    m = Metadata({Metadata::Cell, Metadata::Derived},
+                 std::vector<int>({nhydro + nscalars}), prim_labels);
+  }
   pkg->AddField("prim", m);
 
   const auto refine_str = pin->GetOrAddString("refinement", "type", "unset");

src/hydro/hydro_driver.cpp

@@ -260,6 +260,10 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks,
     auto rkl2_step_first = tl.AddTask(init_MY0, RKL2StepFirst, Y0.get(), base.get(),
                                       Yjm2.get(), MY0.get(), s_rkl, tau);
 
+    // if prolongate on prims is done then this is a ConsToPrim call. Otherwise, noop.
+    auto precomm_fill_derived =
+        tl.AddTask(rkl2_step_first, parthenon::Update::PreCommFillDerived<MeshData<Real>>,
+                   base.get());
     // Update ghost cells of Y1 (as MY1 is calculated for each Y_j).
     // Y1 stored in "base", see rkl2_step_first task.
     // Update ghost cells (local and non local), prolongate and apply bound cond.
@@ -268,7 +272,7 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks,
     // best impl. Go with default call (split local/nonlocal) for now.
     // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars
     auto bounds_exchange = parthenon::AddBoundaryExchangeTasks(
-        rkl2_step_first | start_bnd, tl, base, pmesh->multilevel);
+        precomm_fill_derived | start_bnd, tl, base, pmesh->multilevel);
 
     tl.AddTask(bounds_exchange, parthenon::Update::FillDerived<MeshData<Real>>,
                base.get());
@@ -326,14 +330,18 @@ void AddSTSTasks(TaskCollection *ptask_coll, Mesh *pmesh, BlockList_t &blocks,
           tl.AddTask(set_flx, RKL2StepOther, Y0.get(), base.get(), Yjm2.get(), MY0.get(),
                      mu_j, nu_j, mu_tilde_j, gamma_tilde_j, tau);
 
+      // if prolongate on prims is done then this is a ConsToPrim call. Otherwise, noop.
+      auto precomm_fill_derived =
+          tl.AddTask(rkl2_step_other,
+                     parthenon::Update::PreCommFillDerived<MeshData<Real>>, base.get());
       // update ghost cells of base (currently storing Yj)
       // Update ghost cells (local and non local), prolongate and apply bound cond.
       // TODO(someone) experiment with split (local/nonlocal) comms with respect to
       // performance for various tests (static, amr, block sizes) and then decide on the
       // best impl. Go with default call (split local/nonlocal) for now.
       // TODO(pgrete) optimize (in parthenon) to only send subset of updated vars
       auto bounds_exchange = parthenon::AddBoundaryExchangeTasks(
-          rkl2_step_other | start_bnd, tl, base, pmesh->multilevel);
+          precomm_fill_derived | start_bnd, tl, base, pmesh->multilevel);
 
       tl.AddTask(bounds_exchange, parthenon::Update::FillDerived<MeshData<Real>>,
                  base.get());
@@ -561,11 +569,15 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
           tl.AddTask(source_split_strang_final, AddSplitSourcesFirstOrder, mu0.get(), tm);
     }
 
+    // if prolongate on prims is done then this is a ConsToPrim call. Otherwise, noop.
+    auto precomm_fill_derived =
+        tl.AddTask(source_split_first_order,
+                   parthenon::Update::PreCommFillDerived<MeshData<Real>>, mu0.get());
     // Update ghost cells (local and non local), prolongate and apply bound cond.
     // TODO(someone) experiment with split (local/nonlocal) comms with respect to
     // performance for various tests (static, amr, block sizes) and then decide on the
     // best impl. Go with default call (split local/nonlocal) for now.
-    parthenon::AddBoundaryExchangeTasks(source_split_first_order | start_bnd, tl, mu0,
+    parthenon::AddBoundaryExchangeTasks(precomm_fill_derived | start_bnd, tl, mu0,
                                         pmesh->multilevel);
   }
 

src/pgen/cloud.cpp

@@ -149,9 +149,9 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto hydro_pkg = pmb->packages.Get("Hydro");
-  auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
-  auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
-  auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+  auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire);
+  auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire);
+  auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire);
 
   const auto nhydro = hydro_pkg->Param<int>("nhydro");
   const auto nscalars = hydro_pkg->Param<int>("nscalars");