From f23fc8e10d2df47fc6e5f8bb422f0d40d166f181 Mon Sep 17 00:00:00 2001
From: David Radice <david.radice@psu.edu>
Date: Thu, 8 Aug 2024 16:31:43 +0000
Subject: [PATCH] Fix prolongation with dyn_grmhd

* Ensure that the Z4c variables have been prolongated, before calling
  the ConsToPrim
* Only use Lagrange interpolation to prolongate the Z4c variables
---
 src/driver/driver.cpp        | 33 ++++++++++++++++++---------------
 src/mesh/mesh_refinement.cpp |  7 +++----
 src/mesh/mesh_refinement.hpp |  2 +-
 3 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index fe761ded..0935b380 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -518,6 +518,20 @@ Real Driver::UpdateWallClock() {
 void Driver::InitBoundaryValuesAndPrimitives(Mesh *pm) {
   // Note: with MPI, sends on ALL MBs must be complete before receives execute
 
+  // Initialize Z4c
+  z4c::Z4c *pz4c = pm->pmb_pack->pz4c;
+  if (pz4c != nullptr) {
+    (void) pz4c->RestrictU(this, 0);
+    (void) pz4c->InitRecv(this, -1);  // stage < 0 suppresses InitFluxRecv
+    (void) pz4c->SendU(this, 0);
+    (void) pz4c->ClearSend(this, -1);
+    (void) pz4c->ClearRecv(this, -1);
+    (void) pz4c->RecvU(this, 0);
+    (void) pz4c->Z4cBoundaryRHS(this, 0);
+    (void) pz4c->ApplyPhysicalBCs(this, 0);
+    (void) pz4c->Prolongate(this, 0);
+  }
+
   // Initialize HYDRO: ghost zones and primitive variables (everywhere)
   // includes communications for shearing box boundaries
   hydro::Hydro *phydro = pm->pmb_pack->phydro;
@@ -563,7 +577,10 @@ void Driver::InitBoundaryValuesAndPrimitives(Mesh *pm) {
     if (pdyngr == nullptr) {
       (void) pmhd->ConToPrim(this, 0);
     } else {
-      pdyngr->ConToPrim(this, 0);
+      if (pz4c != nullptr) {
+        (void) pz4c->Z4cToADM_(this, 0);
+      }
+      (void) pdyngr->ConToPrim(this, 0);
     }
   }
 
@@ -581,19 +598,5 @@ void Driver::InitBoundaryValuesAndPrimitives(Mesh *pm) {
     (void) prad->Prolongate(this, 0);
   }
 
-  // Initialize Z4c
-  z4c::Z4c *pz4c = pm->pmb_pack->pz4c;
-  if (pz4c != nullptr) {
-    (void) pz4c->RestrictU(this, 0);
-    (void) pz4c->InitRecv(this, -1);  // stage < 0 suppresses InitFluxRecv
-    (void) pz4c->SendU(this, 0);
-    (void) pz4c->ClearSend(this, -1);
-    (void) pz4c->ClearRecv(this, -1);
-    (void) pz4c->RecvU(this, 0);
-    (void) pz4c->Z4cBoundaryRHS(this, 0);
-    (void) pz4c->ApplyPhysicalBCs(this, 0);
-    (void) pz4c->Prolongate(this, 0);
-  }
-
   return;
 }
diff --git a/src/mesh/mesh_refinement.cpp b/src/mesh/mesh_refinement.cpp
index d9796cca..468640b7 100644
--- a/src/mesh/mesh_refinement.cpp
+++ b/src/mesh/mesh_refinement.cpp
@@ -631,7 +631,7 @@ void MeshRefinement::RedistAndRefineMeshBlocks(ParameterInput *pin, int nnew, in
       RefineFC(new_to_old, pmhd->b0, pmhd->coarse_b0);
     }
     if (pz4c != nullptr) {
-      RefineCC(new_to_old, pz4c->u0, pz4c->coarse_u0);
+      RefineCC(new_to_old, pz4c->u0, pz4c->coarse_u0, true);
     }
   }
 
@@ -1027,11 +1027,10 @@ void MeshRefinement::CopyForRefinementFC(DvceFaceFld4D<Real> &b,DvceFaceFld4D<Re
 //! copied to another location or sent to another rank via MPI.
 
 void MeshRefinement::RefineCC(DualArray1D<int> &n2o, DvceArray5D<Real> &a,
-                              DvceArray5D<Real> &ca) {
+                              DvceArray5D<Real> &ca, bool is_z4c) {
   int nvar = a.extent_int(1);  // TODO(@user): 2nd index from L of in array must be NVAR
   auto &new_nmb = new_nmb_eachrank[global_variable::my_rank];
   MeshBlockPack* pmbp = pmy_mesh->pmb_pack;
-  bool not_z4c = (pmbp->pz4c == nullptr)? true : false;
   auto &indcs = pmy_mesh->mb_indcs;
   auto &cis = indcs.cis, &cie = indcs.cie;
   auto &cjs = indcs.cjs, &cje = indcs.cje;
@@ -1074,7 +1073,7 @@ void MeshRefinement::RefineCC(DualArray1D<int> &n2o, DvceArray5D<Real> &a,
         int fk = 2*k - cks;  // correct when cks=ks
 
         // call inlined prolongation operator for CC variables
-        if (not_z4c) {
+        if (!is_z4c) {
           ProlongCC(m,v,k,j,i,fk,fj,fi,multi_d,three_d,ca,a);
         } else {
           switch (indcs.ng) {
diff --git a/src/mesh/mesh_refinement.hpp b/src/mesh/mesh_refinement.hpp
index 3618f674..2e2f5e6b 100644
--- a/src/mesh/mesh_refinement.hpp
+++ b/src/mesh/mesh_refinement.hpp
@@ -107,7 +107,7 @@ class MeshRefinement {
   void CopyForRefinementCC(DvceArray5D<Real> &a, DvceArray5D<Real> &ca);
   void CopyForRefinementFC(DvceFaceFld4D<Real> &b, DvceFaceFld4D<Real> &cb);
 
-  void RefineCC(DualArray1D<int> &n2o, DvceArray5D<Real> &a, DvceArray5D<Real> &ca);
+  void RefineCC(DualArray1D<int> &n2o, DvceArray5D<Real> &a, DvceArray5D<Real> &ca, bool is_z4c=false);
   void RefineFC(DualArray1D<int> &n2o, DvceFaceFld4D<Real> &b, DvceFaceFld4D<Real> &cb);
 
   void RestrictCC(DvceArray5D<Real> &a, DvceArray5D<Real> &ca, bool is_z4c=false);