Internally consistent boundary logic

Move handling of GRMHD vars at boundaries to Inverter,
try to make comments clearer

kharma/boundaries/boundaries.cpp

@@ -562,7 +562,7 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
     // There are two modes of operation here:
     if (sync_prims) {
         // 1. Exchange/prolongate/restrict PRIMITIVE variables: (ImEx driver)
-        //    Primitive variables and conserved B field are marked FillGhost
+        //    Primitive variables (except B field!) are marked FillGhost
         //    Explicitly run UtoP on B field, then PtoU on everything
         // TODO there should be a set of B field wrappers that dispatch this
         auto pkgs = pmb->packages.AllPackages();

kharma/electrons/electrons.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: electrons.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
@@ -198,6 +198,8 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
 
     pkg->BlockUtoP = Electrons::BlockUtoP;
     pkg->BoundaryUtoP = Electrons::BlockUtoP;
+    // If we wanted to apply the domian boundaries to primitive Electrons variables
+    //pkg->DomainBoundaryPtoU = Electrons::BlockPtoU;
 
     return pkg;
 }
@@ -361,8 +363,8 @@ TaskStatus ApplyElectronHeating(MeshBlockData<Real> *rc_old, MeshBlockData<Real>
             // Heat different electron passives based on different dissipation fraction models
             // Expressions here closely adapted (read: stolen) from implementation in iharm3d
             // courtesy of Cesar Diaz, see https://github.com/AFD-Illinois/iharm3d
-            
-            // In all of these the electron entropy stored value is the entropy conserving solution 
+
+            // In all of these the electron entropy stored value is the entropy conserving solution
                                  // and then when updated it becomes the energy conserving solution
             if (m_p.K_CONSTANT >= 0) {
                 const Real fel = fel_const;
@@ -457,7 +459,7 @@ TaskStatus ApplyElectronHeating(MeshBlockData<Real> *rc_old, MeshBlockData<Real>
             const Real lx2=  pmb->packages.Get("GRMHD")->Param<Real>("lx2");
             const Real edot= pmb->packages.Get("GRMHD")->Param<Real>("drive_edot");
             GridScalar alfven_speed = rc->Get("alfven_speed").data;
-            
+
             int Nx1 = pmb->cellbounds.ncellsi(IndexDomain::interior);
             int Nx2 = pmb->cellbounds.ncellsj(IndexDomain::interior);
             Real *dv0 =  (Real*) malloc(sizeof(Real)*Nx1*Nx2);
@@ -492,9 +494,9 @@ TaskStatus ApplyElectronHeating(MeshBlockData<Real> *rc_old, MeshBlockData<Real>
                     dv0[i*Nx1+j] -= mean_velocity0;
                     dv1[i*Nx1+j] -= mean_velocity1;
                 }
-            } 
+            }
 
-            Real Bhalf = 0; Real A = 0; Real init_e = 0; 
+            Real Bhalf = 0; Real A = 0; Real init_e = 0;
             Kokkos::Sum<Real> Bhalf_reducer(Bhalf); Kokkos::Sum<Real> A_reducer(A); Kokkos::Sum<Real> init_e_reducer(init_e);
             pmb->par_reduce("forced_mhd_normal_kick_normalization_Bhalf", mykb.s, mykb.e, myjb.s, myjb.e, myib.s, myib.e,
                 KOKKOS_LAMBDA(const int k, const int j, const int i, Real &local_result) {
@@ -578,7 +580,7 @@ void ApplyFloors(MeshBlockData<Real> *mbd, IndexDomain domain)
                 } else {
                     ktot_max = floors.ktot_max;
                 }
-                
+
                 if (P(m_p.KTOT, k, j, i) > ktot_max) {
                     fflag(0, k, j, i) = Floors::FFlag::KTOT | (int) fflag(0, k, j, i);
                     P(m_p.KTOT, k, j, i) = ktot_max;

kharma/grmhd/grmhd.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: grmhd.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
@@ -190,16 +190,12 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     // A KHARMAPackage also contains quite a few "callbacks," or functions called at
     // specific points in a step if the package is loaded.
     // Generally, see the headers for function descriptions.
-
-    //pkg->BlockUtoP // Taken care of by separate "Inverter" package since it's hard to do
-
-    // On physical boundaries, even if we've sync'd both, respect the application to primitive variables
-    pkg->DomainBoundaryPtoU = Flux::BlockPtoUMHD;
-
-    // AMR-related
     pkg->CheckRefinementBlock    = GRMHD::CheckRefinement;
     pkg->EstimateTimestepMesh    = GRMHD::EstimateTimestep;
     pkg->PostStepDiagnosticsMesh = GRMHD::PostStepDiagnostics;
+    // NOTE: PtoU and UtoP for the five fluid variables are taken care of by the "Inverter"
+    // package, which registers the relevant callbacks
+    // This is so that the whole package can be disabled when running implicitly or w/EMHD
 
     // List (vector) of HistoryOutputVars that will all be enrolled as output variables
     parthenon::HstVar_list hst_vars = {};

kharma/inverter/inverter.cpp

@@ -107,9 +107,17 @@ std::shared_ptr<KHARMAPackage> Inverter::Initialize(ParameterInput *pin, std::sh
     m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Overridable});
     pkg->AddField("fflag", m);
 
-    // We exist basically to do this
-    pkg->BlockUtoP = Inverter::BlockUtoP;
-    pkg->BoundaryUtoP = Inverter::BlockUtoP;
+    // This package may be loaded even when evolving implicitly, e.g. for FOFC
+    // Only register our callbacks if they're needed for explicit evolution or a guess
+    if (!pin->GetBoolean("GRMHD", "implicit") || pin->GetBoolean("emhd", "ideal_guess")) {
+        // We exist basically to do this
+        pkg->BlockUtoP = Inverter::BlockUtoP;
+        // We want to run U->P on most boundaries when we're synchronizing conserved variables
+        pkg->BoundaryUtoP = Inverter::BlockUtoP;
+        // However, we apply domain boundaries to primitives.
+        // Registering this additional function conveys that to the callers in `Packages` and `Boundaries`
+        pkg->DomainBoundaryPtoU = Flux::BlockPtoUMHD;
+    }
 
     pkg->PostStepDiagnosticsMesh = Inverter::PostStepDiagnostics;