Merge pull request #128 from AFD-Illinois/fix/utop-less

WIP: Stop running unnecessary UtoP at boundaries

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

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: inverter.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
@@ -36,6 +36,7 @@
 // inverter.hpp includes the template and instantiations in the correct order
 
 #include "domain.hpp"
+#include "flux.hpp"
 #include "reductions.hpp"
 
 int Inverter::CountPFlags(MeshData<Real> *md)
@@ -107,9 +108,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;
 
@@ -156,10 +165,10 @@ inline void BlockPerformInversion(MeshBlockData<Real> *rc, IndexDomain domain, b
     const auto& G = pmb->coords;
 
     // Get the primitives from our conserved versions
-    // Notice we recover variables for only the physical (interior or interior-ghost)
+    // Notice by default, we recover variables for only the physical (interior or interior-ghost)
     // zones!  These are the only ones which are filled at our point in the step
-    auto bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
-    const IndexRange3 b = KDomain::GetPhysicalRange(rc);
+    const IndexRange3 b = (domain == IndexDomain::entire)
+                          ? KDomain::GetPhysicalRange(rc) : KDomain::GetRange(rc, domain, coarse);
     pmb->par_for("U_to_P", b.ks, b.ke, b.js, b.je, b.is, b.ie,
         KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
             const Floors::Prescription& myfloors = (inverter_floors.radius_dependent_floors

kharma/kharma_package.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: kharma_package.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
@@ -128,24 +128,24 @@ TaskStatus Packages::BoundaryPtoUElseUtoP(MeshBlockData<Real> *rc, IndexDomain d
     auto pmb = rc->GetBlockPointer();
     auto kpackages = rc->GetBlockPointer()->packages.AllPackagesOfType<KHARMAPackage>();
     // Some downstream UtoP rely on GRMHD prims, some cons
-    if (kpackages.count("GRMHD")) {
-        KHARMAPackage *pkpackage = pmb->packages.Get<KHARMAPackage>("GRMHD");
+    if (kpackages.count("Inverter")) {
+        KHARMAPackage *pkpackage = pmb->packages.Get<KHARMAPackage>("Inverter");
         if (pkpackage->DomainBoundaryPtoU != nullptr) {
-            Flag("DomainBoundaryPtoU_GRMHD");
+            Flag("DomainBoundaryPtoU_Inverter");
             pkpackage->DomainBoundaryPtoU(rc, domain, coarse);
             EndFlag();
         } else if (pkpackage->BoundaryUtoP != nullptr) { // This won't be called
-            Flag("DomainBoundaryUtoP_GRMHD");
+            Flag("DomainBoundaryUtoP_Inverter");
             pkpackage->BoundaryUtoP(rc, domain, coarse);
             EndFlag();
         }
     }
     for (auto kpackage : kpackages) {
-        if (kpackage.second->DomainBoundaryPtoU != nullptr && kpackage.first != "GRMHD") {
+        if (kpackage.second->DomainBoundaryPtoU != nullptr && kpackage.first != "Inverter") {
             Flag("DomainBoundaryPtoU_"+kpackage.first);
             kpackage.second->DomainBoundaryPtoU(rc, domain, coarse);
             EndFlag();
-        } else if (kpackage.second->BoundaryUtoP != nullptr && kpackage.first != "GRMHD") {
+        } else if (kpackage.second->BoundaryUtoP != nullptr && kpackage.first != "Inverter") {
             Flag("DomainBoundaryUtoP_"+kpackage.first);
             kpackage.second->BoundaryUtoP(rc, domain, coarse);
             EndFlag();