Restart from phdf outputs

kharma/b_ct/b_ct.cpp

@@ -94,6 +94,13 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
                                             Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector};
     std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::Conserved, Metadata::WithFluxes,
                                             Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector};
+
+    // KHARMA now (vaguely) supports restarting from dump (.phdf) files.
+    // To do so we need to read cell-centered primitive B and interpolate to faces, as we do with iharm3d restart files
+    if (pin->GetOrAddBoolean("b_field", "restart_from_prims", false)) {
+        flags_prim.push_back(Metadata::Restart);
+    }
+
     std::vector<int> s_vector({NVEC});
     m = Metadata(flags_prim, s_vector);
     pkg->AddField("prims.B", m);
@@ -168,7 +175,17 @@ TaskStatus B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coa
     // Return if we're not syncing U & P at all (e.g. edges)
     if (B_Uf.GetDim(4) == 0) return TaskStatus::complete;
 
-    const IndexRange3 bc = KDomain::GetRange(rc, domain, coarse);
+    // We need one cell inside of the domain, since it will be updated by the domain-face field
+    // this oversteps "interior" by a zone
+    IndexRange3 bct;
+    if (domain == IndexDomain::interior || domain == IndexDomain::entire) {
+        bct = KDomain::GetRange(rc, domain, coarse);
+    } else {
+        const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain);
+        const bool binner = KBoundaries::BoundaryIsInner(bface);
+        bct = KDomain::GetRange(rc, domain, (binner) ? 0 : -1, (binner) ? 1 : 0, coarse);
+    }
+    const IndexRange3 bc = bct;
 
     // Average the primitive vals to zone centers
     pmb->par_for("UtoP_B_center", bc.ks, bc.ke, bc.js, bc.je, bc.is, bc.ie,

kharma/b_flux_ct/b_flux_ct.cpp

@@ -112,6 +112,12 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Independent, Metadata::Restart, Metadata::FillGhost, Metadata::WithFluxes, Metadata::Conserved,
                                             Metadata::Cell, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
 
+    // KHARMA now (vaguely) supports restarting from dump (.phdf) files.
+    // To do so we need to read cell-centered primitive B and interpolate to faces, as we do with iharm3d restart files
+    if (pin->GetOrAddBoolean("b_field", "restart_from_prims", false)) {
+        flags_prim.push_back(Metadata::Restart);
+    }
+
     auto m = Metadata(flags_prim, s_vector);
     pkg->AddField("prims.B", m);
     m = Metadata(flags_cons, s_vector);

kharma/boundaries/boundaries.cpp

@@ -716,6 +716,8 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
                     const Real gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
                     const auto& G = pmb->coords;
                     const int nvar = F.GetDim(4);
+                    // Loci "outer" and "inner" refer to right and left.
+                    // At right edge, our "outer" toward domain is Loci "inner" as in left half
                     const Loci loc = (binner) ? Loci::outer_half : Loci::inner_half;
 
                     const IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, CC);

kharma/driver/kharma_driver.cpp

@@ -1,26 +1,26 @@
 
-/* 
+/*
  *  File: kharma_driver.cpp
- *  
+ *
  *  BSD 3-Clause License
- *  
+ *
  *  Copyright (c) 2020, AFD Group at UIUC
  *  All rights reserved.
- *  
+ *
  *  Redistribution and use in source and binary forms, with or without
  *  modification, are permitted provided that the following conditions are met:
- *  
+ *
  *  1. Redistributions of source code must retain the above copyright notice, this
  *     list of conditions and the following disclaimer.
- *  
+ *
  *  2. Redistributions in binary form must reproduce the above copyright notice,
  *     this list of conditions and the following disclaimer in the documentation
  *     and/or other materials provided with the distribution.
- *  
+ *
  *  3. Neither the name of the copyright holder nor the names of its
  *     contributors may be used to endorse or promote products derived from
  *     this software without specific prior written permission.
- *  
+ *
  *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -85,8 +85,8 @@ std::shared_ptr<KHARMAPackage> KHARMADriver::Initialize(ParameterInput *pin, std
     // Add a flag if we wish to use ideal variables explicitly evolved as guess for implicit update.
     // GRIM uses the fluid state of the previous (sub-)step.
     // The logic here is that the non-ideal variables do not significantly contribute to the
-    // stress-energy tensor. We can explicitly update the ideal MHD variables and hope that this 
-    // takes us to a region in the parameter space of state variables that 
+    // stress-energy tensor. We can explicitly update the ideal MHD variables and hope that this
+    // takes us to a region in the parameter space of state variables that
     // is close to the true solution. The corrections obtained from the implicit update would then
     // be small corrections.
     Metadata::AddUserFlag("IdealGuess");
@@ -187,14 +187,14 @@ TaskStatus KHARMADriver::SyncAllBounds(std::shared_ptr<MeshData<Real>> &md)
     Flag("SyncAllBounds");
     TaskID t_none(0);
 
-    //MPIBarrier();
+    MPIBarrier();
 
     TaskCollection tc;
     auto tr = tc.AddRegion(1);
     AddBoundarySync(t_none, tr[0], md);
     while (!tr.Execute());
 
-    //MPIBarrier();
+    MPIBarrier();
 
     EndFlag();
     return TaskStatus::complete;

kharma/prob/post_initialize.cpp

@@ -88,6 +88,8 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
             // normalize the magnetic field according to the max density
             bool is_torus = prob_name == "torus";
             if (pin->GetOrAddBoolean("b_field", "norm", is_torus)) {
+                // Sync to apply physical boundary conditions and fill B ghosts correctly
+                KHARMADriver::SyncAllBounds(md);
                 NormalizeBField(md.get(), pin);
             }
         }
@@ -112,7 +114,10 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         // We only record the conserved magnetic field in KHARMA restarts,
         // but we record primitive field in iharm3d restarts
         bool iharm3d_restart = prob_name == "resize_restart";
-        if (!iharm3d_restart) {
+        // but sometimes in a pinch we want to restart from .phdf files, which are like iharm3d restarts
+        bool prims_only_restart = pin->GetBoolean("b_field", "restart_from_prims");
+        if (!(iharm3d_restart  || prims_only_restart)) {
+            std::cerr << "Restoring B field from conserved values" << std::endl;
             if (pkgs.count("B_FluxCT")) {
                 B_FluxCT::MeshUtoP(md.get(), IndexDomain::entire);
             } else if (pkgs.count("B_CT")) {
@@ -122,6 +127,7 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
                 EMHD::MeshUtoP(md.get(), IndexDomain::entire);
             }
         } else {
+            std::cerr << "Restoring B field from primitive values" << std::endl;
             if (pkgs.count("B_FluxCT")) {
                 B_FluxCT::MeshPtoU(md.get(), IndexDomain::entire);
             } else if (pkgs.count("B_CT")) {

kharma/prob/seed_B.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: seed_B.cpp
- *  
+ *
  *  BSD 3-Clause License
- *  
+ *
  *  Copyright (c) 2020, AFD Group at UIUC
  *  All rights reserved.
- *  
+ *
  *  Redistribution and use in source and binary forms, with or without
  *  modification, are permitted provided that the following conditions are met:
- *  
+ *
  *  1. Redistributions of source code must retain the above copyright notice, this
  *     list of conditions and the following disclaimer.
- *  
+ *
  *  2. Redistributions in binary form must reproduce the above copyright notice,
  *     this list of conditions and the following disclaimer in the documentation
  *     and/or other materials provided with the distribution.
- *  
+ *
  *  3. Neither the name of the copyright holder nor the names of its
  *     contributors may be used to endorse or promote products derived from
  *     this software without specific prior written permission.
- *  
+ *
  *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -47,19 +47,20 @@ using namespace parthenon;
 
 /**
  * Perform a Parthenon MPI reduction.
- * Should only be used in initialization code, as the
- * reducer object & MPI comm are created on entry &
- * cleaned on exit
+ * Should only be used in initialization code, as it
+ * is very not fast.
  */
 template<typename T>
 inline T MPIReduce_once(T f, MPI_Op O)
 {
+    MPIBarrier();
     parthenon::AllReduce<T> reduction;
     reduction.val = f;
     reduction.StartReduce(O);
     // Wait on results
     while (reduction.CheckReduce() == parthenon::TaskStatus::incomplete);
     // TODO catch errors?
+    MPIBarrier();
     return reduction.val;
 }
 
@@ -114,7 +115,7 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
     if constexpr (Seed == BSeedType::constant ||
                   Seed == BSeedType::monopole ||
                   Seed == BSeedType::orszag_tang ||
-                  Seed == BSeedType::wave || 
+                  Seed == BSeedType::wave ||
                   Seed == BSeedType::shock_tube)
     {
         // All custom B fields should set what they need of these.
@@ -331,7 +332,8 @@ TaskStatus SeedBFieldType(MeshBlockData<Real> *rc, ParameterInput *pin, IndexDom
                     // So, we preserve exact values in the no-tilt case.
                     A(V3, k, j, i) = Aphi;
                 }
-            });
+            }
+        );
 
         if (pkgs.count("B_CT")) {
             auto B_Uf = rc->PackVariables(std::vector<std::string>{"cons.fB"});
@@ -526,9 +528,6 @@ TaskStatus NormalizeBField(MeshData<Real> *md, ParameterInput *pin)
         }
     }
 
-    // We've been initializing/manipulating P
-    Flux::MeshPtoU(md, IndexDomain::entire);
-
     EndFlag(); //NormBField
     return TaskStatus::complete;
 }