Merge pull request #132 from AFD-Illinois/fix/papercuts

Polar Papercuts

kharma/b_cd/b_cd.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: b_cd.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
@@ -246,7 +246,7 @@ TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData<Real> *md)
 
 void FillOutput(MeshBlock *pmb, ParameterInput *pin)
 {
-    auto rc = pmb->meshblock_data.Get().get();
+    auto rc = pmb->meshblock_data.Get("base").get();
     IndexDomain domain = IndexDomain::interior;
     int is = pmb->cellbounds.is(domain), ie = pmb->cellbounds.ie(domain);
     int js = pmb->cellbounds.js(domain), je = pmb->cellbounds.je(domain);

kharma/b_ct/b_ct.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: b_ct.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
@@ -582,7 +582,8 @@ void B_CT::CalcDivB(MeshData<Real> *md, std::string divb_field_name)
 
 void B_CT::FillOutput(MeshBlock *pmb, ParameterInput *pin)
 {
-    auto rc = pmb->meshblock_data.Get();
+    // This is called after the step, use "base" container
+    auto rc = pmb->meshblock_data.Get("base");
     const int ndim = pmb->pmy_mesh->ndim;
     if (ndim < 2) return;
 

kharma/b_ct/b_ct_boundaries.cpp

@@ -35,8 +35,10 @@
 
 #include "decs.hpp"
 #include "domain.hpp"
+#include "floors.hpp"
 #include "grmhd.hpp"
 #include "grmhd_functions.hpp"
+#include "inverter.hpp"
 #include "kharma.hpp"
 
 void B_CT::ZeroBoundaryEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &emfpack, bool coarse)
@@ -285,12 +287,21 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
     const int bdir = KBoundaries::BoundaryDirection(bface);
     const auto bname = KBoundaries::BoundaryName(bface);
 
+    // Pull cell-centered values, as we need to update fluid primitives
     // TODO standardize on passing Packs or Datas...
-    auto B_U = rc->PackVariables(std::vector<std::string>{"cons.B"});
-    auto B_P = rc->PackVariables(std::vector<std::string>{"prims.B"});
+    PackIndexMap prims_map, cons_map;
+    auto P = rc->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map);
+    auto U = rc->PackVariables(std::vector<MetadataFlag>{Metadata::Conserved, Metadata::Cell}, cons_map);
+    const VarMap m_u(cons_map, true), m_p(prims_map, false);
 
     const auto& G = pmb->coords;
 
+    const Real gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
+
+    const Floors::Prescription floors = pmb->packages.Get("Floors")->Param<Floors::Prescription>("prescription");
+    // Don't be fooled, this function does *not* support/preserve EMHD values
+    const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb->packages);
+
     // Subtract the average B3 as "reconnection"
     IndexRange3 b = KDomain::GetRange(rc, domain, F3, coarse);
     IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, F3, coarse);
@@ -318,12 +329,22 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
             );
             member.team_barrier();
 
-            // Update cell-centered conserved & primitive B3. Not worth a separate BlockUtoP call
-            parthenon::par_for_inner(member, b.ks, b.ke-1,
+            // Update cell-centered conserved & primitive B, and cell primitive fluid variables, in the zones we touched
+            parthenon::par_for_inner(member, b.ks, b.ke-1, // iterate over all *cells* k
                 [&](const int& k) {
-                    B_P(V3, k, jf, i) =  (fpack(F3, 0, k, jf, i) / G.gdet(Loci::face3, jf, i)
+                    P(m_p.B3, k, jf, i) =  (fpack(F3, 0, k, jf, i) / G.gdet(Loci::face3, jf, i)
                                         + fpack(F3, 0, k + 1, jf, i) / G.gdet(Loci::face3, jf, i)) / 2;
-                    B_U(V3, k, jf, i) = B_P(V3, k, jf, i) * G.gdet(Loci::center, jf, i);
+                    U(m_u.B3, k, jf, i) = P(m_p.B3, k, jf, i) * G.gdet(Loci::center, jf, i);
+
+                    // Recover primitive GRMHD variables from our modified U
+                    Inverter::u_to_p<Inverter::Type::kastaun>(G, U, m_u, gam, k, jf, i, P, m_p, Loci::center,
+                                                              floors, 8, 1e-8);
+                    // Floor them
+                    int fflag = Floors::apply_geo_floors(G, P, m_p, gam, k, jf, i, floors, floors, Loci::center);
+                    // Recalculate U on anything we floored
+                    if (fflag)
+                        GRMHD::p_to_u(G, P, m_p, gam, k, jf, i, U, m_u, Loci::center);
+
                 }
             );
         }

kharma/b_flux_ct/b_flux_ct.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: b_flux_ct.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
@@ -410,7 +410,7 @@ void ZeroBoundaryFlux(MeshData<Real> *md, IndexDomain domain, bool coarse)
     // Compare this section with calculation of emf3 in FluxCT:
     // these changes ensure that boundary emfs emf3(i,js,k)=0, etc.
     for (auto &pmb : pmesh->block_list) {
-        auto& rc = pmb->meshblock_data.Get();
+        auto& rc = pmb->meshblock_data.Get(md->StageName());
         auto& B_F = rc->PackVariablesAndFluxes(std::vector<std::string>{"cons.B"});
 
         if (domain == IndexDomain::inner_x2 &&
@@ -495,7 +495,7 @@ void Bflux0(MeshData<Real> *md, IndexDomain domain, bool coarse)
     // Compare this section with calculation of emf3 in FluxCT:
     // these changes ensure that boundary emfs emf3(i,js,k)=0, etc.
     for (auto &pmb : pmesh->block_list) {
-        auto& rc = pmb->meshblock_data.Get();
+        auto& rc = pmb->meshblock_data.Get(md->StageName());
         auto& B_F = rc->PackVariablesAndFluxes(std::vector<std::string>{"cons.B"});
 
         // "Bflux0" prescription for keeping divB~=0 on zone corners of the interior & exterior X1 faces
@@ -555,7 +555,7 @@ IndexRange ValidDivBX1(MeshBlock *pmb)
 {
     // All user, physical (not MPI/periodic) boundary conditions in X1 will generate divB on corners
     // intersecting the interior & exterior faces. Don't report these zones, as we expect it.
-    const IndexRange ibl = pmb->meshblock_data.Get()->GetBoundsI(IndexDomain::interior);
+    const IndexRange ibl = pmb->meshblock_data.Get("base")->GetBoundsI(IndexDomain::interior);
     bool avoid_inner = (!pmb->packages.Get("B_FluxCT")->Param<bool>("fix_flux_inner_x1") &&
         pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user);
     bool avoid_outer = (!pmb->packages.Get("B_FluxCT")->Param<bool>("fix_flux_outer_x1") &&
@@ -677,7 +677,7 @@ void CalcDivB(MeshData<Real> *md, std::string divb_field_name)
 }
 void FillOutput(MeshBlock *pmb, ParameterInput *pin)
 {
-    auto rc = pmb->meshblock_data.Get().get();
+    auto rc = pmb->meshblock_data.Get("base").get();
     const int ndim = pmb->pmy_mesh->ndim;
     if (ndim < 2) return;
 

kharma/boundaries/boundaries.cpp

@@ -192,7 +192,7 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
             // introduce divergence to the first physical zone.
             bool clean_face_B = pin->GetOrAddBoolean("boundaries", "clean_face_B_" + bname, (btype == "outflow"));
             params.Add("clean_face_B_"+bname, clean_face_B);
-            // Forcibly reconnect field loops that get trapped around the polar boundary.  Probably not needed anymore.
+            // Forcibly reconnect field loops that get trapped around the polar boundary
             bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, false);
             params.Add("reconnect_B3_"+bname, reconnect_B3);
 
@@ -626,7 +626,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
     const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, -1, 1);
 
     for (auto &pmb : pmesh->block_list) {
-        auto &rc = pmb->meshblock_data.Get();
+        auto &rc = pmb->meshblock_data.Get(md->StageName());
 
         for (int i = 0; i < BOUNDARY_NFACES; i++) {
             BoundaryFace bface = (BoundaryFace)i;

kharma/coord_output/coord_output.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: coord_output.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
@@ -70,7 +70,7 @@ std::shared_ptr<KHARMAPackage> CoordinateOutput::Initialize(ParameterInput *pin,
     pkg->AddField("coords.r", m0);
     pkg->AddField("coords.th", m0);
     pkg->AddField("coords.phi", m0);
-    
+
     // Metric
     pkg->AddField("coords.gcon", m2);
     pkg->AddField("coords.gcov", m2);
@@ -100,7 +100,7 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter
 {
     auto& globals = pmb->packages.Get("Globals")->AllParams();
     if (!globals.Get<bool>("in_loop")) {
-        auto rc = pmb->meshblock_data.Get();
+        auto rc = pmb->meshblock_data.Get("base");
 
         PackIndexMap geom_map;
         auto Geom = rc->PackVariables({Metadata::GetUserFlag("CoordinateOutput")}, geom_map);

kharma/current/current.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: current.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
@@ -118,7 +118,7 @@ void Current::FillOutput(MeshBlock *pmb, ParameterInput *pin)
 {
     // The "preserve" container will only exist after we've taken a step,
     // catch that situation
-    auto& rc1 = pmb->meshblock_data.Get();
+    auto& rc1 = pmb->meshblock_data.Get("base");
     auto rc0 = rc1; // Avoid writing rc0's type when initializing. Still light.
     try {
         // Get the state at beginning of the step

kharma/driver/simple_step.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: simple_step.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
@@ -52,7 +52,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s
 
     // Allocate the fluid states ("containers") we need for each block
     for (auto& pmb : blocks) {
-        auto &base = pmb->meshblock_data.Get();
+        auto &base = pmb->meshblock_data.Get("base");
         if (stage == 1) {
             pmb->meshblock_data.Add("dUdt", base);
             for (int i = 1; i < integrator->nstages; i++)
@@ -68,7 +68,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s
     TaskRegion &single_tasklist_per_pack_region = tc.AddRegion(num_partitions);
     for (int i = 0; i < num_partitions; i++) {
         auto &tl = single_tasklist_per_pack_region[i];
-        // Container names: 
+        // Container names:
         // '_full_step_init' refers to the fluid state at the start of the full time step (Si in iharm3d)
         // '_sub_step_init' refers to the fluid state at the start of the sub step (Ss in iharm3d)
         // '_sub_step_final' refers to the fluid state at the end of the sub step (Sf in iharm3d)
@@ -86,7 +86,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s
         // Any package modifications to the fluxes.  e.g.:
         // 1. CT calculations for B field transport
         // 2. Zero fluxes through poles
-        // etc 
+        // etc
         auto t_fix_flux = tl.AddTask(t_fluxes, Packages::FixFlux, md_sub_step_init.get());
 
         // Apply the fluxes to calculate a change in cell-centered values "md_flux_src"