Expand some indices when checking for inflows, make GetRange usable f…

…or boundary ranges
AFD-Illinois · Nov 13, 2023 · ece0e53 · ece0e53
1 parent 546594e
commit ece0e53
Show file tree

Hide file tree

Showing 6 changed files with 104 additions and 18 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,7 @@
+# Memory dumps
+core*
+*.swp
+
 # Various script results/logs
 out-*.txt
 *.json
@@ -6,7 +10,6 @@ convergence.txt
 *.png
 *.mp4
 *.webm
-core.*
 frames_*/
 logs/
 *.log
@@ -77,6 +80,3 @@ make_args
 # Python files
 __pycache__/
 *.pyc
-
-# added by Hyerin
-*.swp
diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp
@@ -378,9 +378,20 @@ void KBoundaries::CheckInflow(std::shared_ptr<MeshBlockData<Real>> &rc, IndexDom
 
     // Inflow check
     // Iterate over zones w/p=0
-    pmb->par_for_bndry(
-        "check_inflow", IndexRange{0, 0}, domain, CC, coarse,
-        KOKKOS_LAMBDA(const int &p, const int &k, const int &j, const int &i) {
+    // pmb->par_for_bndry(
+    //     "check_inflow", IndexRange{0, 0}, domain, CC, coarse,
+    //     KOKKOS_LAMBDA(const int &p, const int &k, const int &j, const int &i) {
+    //         KBoundaries::check_inflow(G, P, domain, m_p.U1, k, j, i);
+    //     }
+    // );
+    const auto bface = BoundaryFace(domain);
+    const auto bname = BoundaryName(bface);
+    const bool binner = BoundaryIsInner(bface);
+    // One domain interior to boundary
+    auto b = KDomain::GetRange(rc, domain, -((int) !binner), binner, coarse);
+    pmb->par_for(
+        "zero_inflow_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie,
+        KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
             KBoundaries::check_inflow(G, P, domain, m_p.U1, k, j, i);
         }
     );
@@ -427,6 +438,10 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
     // These functions do *not* need an extra row outside the domain,
     // like B_FluxCT::FixBoundaryFlux does.
     const int ndim = pmesh->ndim;
+    // Entire range
+    const IndexRange ibe = pmb0->cellbounds.GetBoundsI(IndexDomain::entire);
+    const IndexRange jbe = pmb0->cellbounds.GetBoundsJ(IndexDomain::entire);
+    const IndexRange kbe = pmb0->cellbounds.GetBoundsK(IndexDomain::entire);
     // Ranges for sides
     const IndexRange ibs = pmb0->cellbounds.GetBoundsI(IndexDomain::interior);
     const IndexRange jbs = pmb0->cellbounds.GetBoundsJ(IndexDomain::interior);
@@ -448,7 +463,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
             if (bdir > ndim) continue;
 
             // Set ranges based
-            IndexRange ib = ibs, jb = jbs, kb = kbs;
+            IndexRange ib = ibe, jb = jbe, kb = kbe;
             // Range for inner_x1 bounds is first face only, etc.
             if (bdir == 1) {
                 ib.s = ib.e = (binner) ? ibf.s : ibf.e;
@@ -468,14 +483,14 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
                 if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
                     if (binner) {
                         pmb->par_for(
-                            "zero_inflow_flux_" + bname, kb.s, kb.e, jb.s, jb.e, ib.s, ib.s,
+                            "zero_inflow_flux_" + bname, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
                             KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
                                 F.flux(bdir, m_rho, k, j, i) = m::min(F.flux(bdir, m_rho, k, j, i), 0.);
                             }
                         );
                     } else {
                         pmb->par_for(
-                            "zero_inflow_flux_" + bname, kb.s, kb.e, jb.s, jb.e, ib.s, ib.s,
+                            "zero_inflow_flux_" + bname, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
                             KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
                                 F.flux(bdir, m_rho, k, j, i) = m::max(F.flux(bdir, m_rho, k, j, i), 0.);
                             }

diff --git a/kharma/domain.hpp b/kharma/domain.hpp
@@ -110,14 +110,18 @@ inline IndexRange3 GetRange(T data, IndexDomain domain, int left_halo=0, int rig
     const IndexRange jb = cellbounds.GetBoundsJ(domain);
     const IndexRange kb = cellbounds.GetBoundsK(domain);
     // Compute sizes with specified halo zones included in non-trivial dimensions
+    // TODO notion of activated x1+x3 with nx2==0?
     const int& ndim = GetNDim(data);
-    // If ghost & not x1 direction
     const IndexRange il = IndexRange{ib.s + left_halo, ib.e + right_halo};
     const IndexRange jl = (ndim > 1) ? IndexRange{jb.s + left_halo, jb.e + right_halo} : jb;
     const IndexRange kl = (ndim > 2) ? IndexRange{kb.s + left_halo, kb.e + right_halo} : kb;
-    return IndexRange3{(uint) il.s, (uint) il.e,
-                       (uint) jl.s, (uint) jl.e,
-                       (uint) kl.s, (uint) kl.e};
+    // Bounds of entire domain, we never mean to go beyond these
+    const IndexRange ibe = cellbounds.GetBoundsI(IndexDomain::entire);
+    const IndexRange jbe = cellbounds.GetBoundsJ(IndexDomain::entire);
+    const IndexRange kbe = cellbounds.GetBoundsK(IndexDomain::entire);
+    return IndexRange3{(uint) m::max(il.s, ibe.s), (uint) m::min(il.e, ibe.e),
+                       (uint) m::max(jl.s, jbe.s), (uint) m::min(jl.e, jbe.e),
+                       (uint) m::max(kl.s, kbe.s), (uint) m::min(kl.e, kbe.e)};
 }
 template<typename T>
 inline IndexRange3 GetRange(T data, IndexDomain domain, bool coarse)

diff --git a/kharma/floors/floors.cpp b/kharma/floors/floors.cpp
@@ -166,8 +166,8 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData<Real> *
     auto pmb = mbd->GetBlockPointer();
 
     PackIndexMap prims_map, cons_map;
-    auto P = mbd->PackVariables({Metadata::GetUserFlag("Primitive")}, prims_map);
-    auto U = mbd->PackVariables(std::vector<MetadataFlag>{Metadata::Conserved}, cons_map);
+    auto P = mbd->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map);
+    auto U = mbd->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;

diff --git a/pars/benchmark/sane_perf_emhd.par b/pars/benchmark/sane_perf_emhd.par
@@ -0,0 +1,67 @@
+# SANE model emulating a real run for performance testing
+# Takes only 1k steps, dramatically reduced dump files
+# Uses HARM driver harm_driver.cpp
+
+# (Also no archival parfile, B cleanup, or two-sync)
+
+<parthenon/job>
+problem_id = torus
+
+# 8 meshblocks -> up to 2 nodes.
+# Pretty representative size for a long simulation
+# Larger simulations have smaller in-simulation timesteps
+<parthenon/mesh>
+refinement = none
+numlevel = 1
+nx1 = 256
+nx2 = 128
+nx3 = 256
+
+<parthenon/meshblock>
+nx1 = 128
+nx2 = 32
+nx3 = 128
+
+<coordinates>
+base = spherical_ks
+transform = fmks
+r_out = 1000
+a = 0.9375
+
+<parthenon/time>
+tlim = 10000.0
+# Limit to 1k steps
+nlim = 1000
+
+<GRMHD>
+cfl = 0.8
+gamma = 1.666667
+reconstruction = weno5
+
+
+
+<driver>
+type = kharma
+two_sync = true
+
+<torus>
+rin = 6.0
+rmax = 12.0
+
+<perturbation>
+u_jitter = 0.04
+
+<b_field>
+type = sane
+beta_min = 100.
+
+<floors>
+rho_min_geom = 1e-6
+u_min_geom = 1e-8
+bsq_over_rho_max = 100
+u_over_rho_max = 2
+
+<debug>
+verbose = 1
+extra_checks = 1
+flag_verbose = 0
diff --git a/pars/tori_3d/sane.par b/pars/tori_3d/sane.par
@@ -20,7 +20,7 @@ nx3 = 32
 <coordinates>
 base = spherical_ks
 transform = fmks
-r_out = 1000
+r_out = 100
 a = 0.9375
 hslope = 0.3
 mks_smooth = 0.5
@@ -71,7 +71,7 @@ Tp = 10
 file_type = hdf5
 dt = 5.0
 single_precision_output = true
-variables = prims.rho, prims.u, prims.uvec, prims.B, jcon, divB
+variables = prims, jcon, divB
 
 <parthenon/output1>
 file_type = rst