diff --git a/.gitignore b/.gitignore index 2e87022a..0cf1683b 100644 --- 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 index a7afe62e..8d08c468 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -378,9 +378,20 @@ void KBoundaries::CheckInflow(std::shared_ptr> &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 *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 *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 *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 index c100b3a2..b3cc6ee3 100644 --- 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 inline IndexRange3 GetRange(T data, IndexDomain domain, bool coarse) diff --git a/kharma/floors/floors.cpp b/kharma/floors/floors.cpp index e1a7ee76..e92aff7d 100644 --- a/kharma/floors/floors.cpp +++ b/kharma/floors/floors.cpp @@ -166,8 +166,8 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData * auto pmb = mbd->GetBlockPointer(); PackIndexMap prims_map, cons_map; - auto P = mbd->PackVariables({Metadata::GetUserFlag("Primitive")}, prims_map); - auto U = mbd->PackVariables(std::vector{Metadata::Conserved}, cons_map); + auto P = mbd->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map); + auto U = mbd->PackVariables(std::vector{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 new file mode 100644 index 00000000..94f06a86 --- /dev/null +++ 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) + + +problem_id = torus + +# 8 meshblocks -> up to 2 nodes. +# Pretty representative size for a long simulation +# Larger simulations have smaller in-simulation timesteps + +refinement = none +numlevel = 1 +nx1 = 256 +nx2 = 128 +nx3 = 256 + + +nx1 = 128 +nx2 = 32 +nx3 = 128 + + +base = spherical_ks +transform = fmks +r_out = 1000 +a = 0.9375 + + +tlim = 10000.0 +# Limit to 1k steps +nlim = 1000 + + +cfl = 0.8 +gamma = 1.666667 +reconstruction = weno5 + + + + +type = kharma +two_sync = true + + +rin = 6.0 +rmax = 12.0 + + +u_jitter = 0.04 + + +type = sane +beta_min = 100. + + +rho_min_geom = 1e-6 +u_min_geom = 1e-8 +bsq_over_rho_max = 100 +u_over_rho_max = 2 + + +verbose = 1 +extra_checks = 1 +flag_verbose = 0 diff --git a/pars/tori_3d/sane.par b/pars/tori_3d/sane.par index 08d94b4e..2d2656ab 100644 --- a/pars/tori_3d/sane.par +++ b/pars/tori_3d/sane.par @@ -20,7 +20,7 @@ nx3 = 32 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 file_type = rst