Merge pull request #575 from IAS-Astrophysics/dyngrmhd-refine-fix

Fix prolongation with dyn_grmhd
IAS-Astrophysics · Aug 16, 2024 · 826fbd9 · 826fbd9
2 parents 3df8cbc + 75b51ef
commit 826fbd9
Show file tree

Hide file tree

Showing 4 changed files with 25 additions and 23 deletions.
diff --git 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
@@ -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) {
@@ -1184,7 +1183,6 @@ void MeshRefinement::RestrictCC(DvceArray5D<Real> &u, DvceArray5D<Real> &cu,
   int nvar = u.extent_int(1);  // TODO(@user): 2nd index from L of in array must be NVAR
 
   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;
@@ -1219,7 +1217,7 @@ void MeshRefinement::RestrictCC(DvceArray5D<Real> &u, DvceArray5D<Real> &cu,
       int finei = 2*i - cis;  // correct when cis=is
       int finej = 2*j - cjs;  // correct when cjs=js
       int finek = 2*k - cks;  // correct when cks=ks
-      if (not_z4c) {
+      if (!is_z4c) {
         cu(m,n,k,j,i) =
             0.125*(u(m,n,finek  ,finej  ,finei) + u(m,n,finek  ,finej  ,finei+1)
                 + u(m,n,finek  ,finej+1,finei) + u(m,n,finek  ,finej+1,finei+1)

diff --git a/src/mesh/mesh_refinement.hpp b/src/mesh/mesh_refinement.hpp
@@ -107,7 +107,8 @@ 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);

diff --git a/src/z4c/z4c_tasks.cpp b/src/z4c/z4c_tasks.cpp
@@ -421,7 +421,7 @@ TaskStatus Z4c::RestrictWeyl(Driver *pdrive, int stage) {
     float time_32 = static_cast<float>(pmy_pack->pmesh->time);
     if ((last_output_time==time_32) && (stage == pdrive->nexp_stages)) {
       if (pmy_pack->pmesh->multilevel) {
-        pmy_pack->pmesh->pmr->RestrictCC(u_weyl, coarse_u_weyl);
+        pmy_pack->pmesh->pmr->RestrictCC(u_weyl, coarse_u_weyl, true);
       }
     }
     return TaskStatus::complete;