diff --git a/external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp b/external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp index 9552fa06..e7374341 100644 --- a/external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp +++ b/external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp @@ -5,8 +5,6 @@ #include "KokkosBatched_Util.hpp" -#define KOKKOS_IMPL_DO_NOT_USE_PRINTF(...) ::printf(__VA_ARGS__) - namespace KokkosBatched { /// diff --git a/kharma/b_cd/b_cd.cpp b/kharma/b_cd/b_cd.cpp index 6e917066..f772dcdd 100644 --- a/kharma/b_cd/b_cd.cpp +++ b/kharma/b_cd/b_cd.cpp @@ -38,13 +38,15 @@ using namespace parthenon; +// THIS MODULE DOESN'T WORK with KHARMA's initialization or modern structure +// The code is here so that we can ensure it keeps compiling, +// which should make it easier to reintroduce if we want to later + namespace B_CD { std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr& packages) { - throw std::runtime_error("Constraint-damping transport is not functional with modern B field initialization!"); - auto pkg = std::make_shared("B_CD"); Params ¶ms = pkg->AllParams(); @@ -99,7 +101,12 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

AddParam<>(parthenon::hist_param_key, hst_vars); - return pkg; + // Throw down here, like this, to avoid inaccessible code warnings + if (1) { + throw std::runtime_error("Constraint-damping transport is not functional with modern B field initialization!"); + } else { + return pkg; + } } void BlockUtoP(MeshBlockData *rc, IndexDomain domain, bool coarse) diff --git a/kharma/b_cd/b_cd.hpp b/kharma/b_cd/b_cd.hpp index eafd3820..c0c4ebf8 100644 --- a/kharma/b_cd/b_cd.hpp +++ b/kharma/b_cd/b_cd.hpp @@ -43,6 +43,10 @@ using namespace parthenon; /** + * THIS MODULE DOESN'T WORK with KHARMA's initialization or modern structure + * The code is here so that we can ensure it keeps compiling, + * which should make it easier to reintroduce if we want to later + * * This physics package implements B field transport with Constraint-Damping (Dedner et al 2002) * * This requires only the values at cell centers, and preserves a cell-centered divergence representation diff --git a/kharma/b_cleanup/b_cleanup.cpp b/kharma/b_cleanup/b_cleanup.cpp index 51276f26..54594234 100644 --- a/kharma/b_cleanup/b_cleanup.cpp +++ b/kharma/b_cleanup/b_cleanup.cpp @@ -139,11 +139,13 @@ std::shared_ptr B_Cleanup::Initialize(ParameterInput *pin, std::s int cleanup_interval = pin->GetOrAddInteger("b_cleanup", "cleanup_interval", manage_field ? 10 : -1); params.Add("cleanup_interval", cleanup_interval); - // Declare fields if we're doing that if (manage_field) { - // Stolen verbatim from FluxCT, will need updates to actually use - throw std::runtime_error("B field cleanup/projection is set as B field transport! If you really want this, disable this error in source!"); + throw std::runtime_error("B field cleanup/projection is set as B field transport! If you really want this, disable this error in source code!"); + } + // Declare fields if we're doing that + if (manage_field) { + // Stolen verbatim from FluxCT, will need updates from there to actually use // Mark if we're evolving implicitly bool implicit_b = pin->GetOrAddBoolean("b_field", "implicit", false); params.Add("implicit", implicit_b); diff --git a/kharma/b_ct/b_ct.cpp b/kharma/b_ct/b_ct.cpp index 864ace84..086341d7 100644 --- a/kharma/b_ct/b_ct.cpp +++ b/kharma/b_ct/b_ct.cpp @@ -38,8 +38,6 @@ #include "grmhd.hpp" #include "grmhd_functions.hpp" #include "kharma.hpp" -// TODO eliminate sync -#include "kharma_driver.hpp" #include #include @@ -71,7 +69,7 @@ std::shared_ptr B_CT::Initialize(ParameterInput *pin, std::shared // TODO gs05_alpha, LDZ04 UCT1, LDZ07 UCT2 std::vector ct_scheme_options = {"bs99", "gs05_0", "gs05_c", "sg07"}; - std::string ct_scheme = pin->GetOrAddString("b_field", "ct_scheme", "bs99", ct_scheme_options); + std::string ct_scheme = pin->GetOrAddString("b_field", "ct_scheme", "sg07", ct_scheme_options); params.Add("ct_scheme", ct_scheme); if (ct_scheme == "gs05_c") std::cout << "KHARMA WARNING: G&S '05 epsilon_c CT is not well-tested." << std::endl @@ -445,132 +443,6 @@ TaskStatus B_CT::AddSource(MeshData *md, MeshData *mdudt, IndexDomai return TaskStatus::complete; } -void B_CT::ZeroEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse) -{ - // TODO might be able to get away with not zeroing ghost EMFs, - // since they shouldn't be used. - auto pmb = rc->GetBlockPointer(); - const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); - const std::string bname = KBoundaries::BoundaryName(bface); - for (auto &el : {E1, E2, E3}) { - // This inlcudes e.g. E1/E3 edges on X2 face - auto b = KDomain::GetBoundaryRange(rc, domain, el, coarse); - pmb->par_for( - "zero_EMF_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - emfpack(el, 0, k, j, i) = 0; - } - ); - } -} - -void B_CT::AverageEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse) -{ - auto pmb = rc->GetBlockPointer(); - const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); - const std::string bname = KBoundaries::BoundaryName(bface); - const int bdir = KBoundaries::BoundaryDirection(bface); - const bool binner = KBoundaries::BoundaryIsInner(bface); - if (bdir != 2) { - throw std::runtime_error("Polar average EMF implemented only in X2!"); - } - - // X1 and X2 EMF are zeroed only *within* boundary domain - // TODO might be able to get away with not zeroing ghost EMFs, - // since they shouldn't be used. - for (auto el : {E1, E2}) { - IndexRange3 b = KDomain::GetRange(rc, domain, el, coarse); - pmb->par_for( - "zero_offaxis_EMF12_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - emfpack(el, 0, k, j, i) = 0.; - } - ); - } - // X3 EMF must additionally be zero *on* polar face, since edge size is 0 - IndexRange3 b = KDomain::GetBoundaryRange(rc, domain, E3, coarse); - pmb->par_for( - "zero_offaxis_EMF3_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - emfpack(E3, 0, k, j, i) = 0; - } - ); - - // In 2D, "averaging" should just mean not zeroing E1 - if (KDomain::GetNDim(rc) < 3) return; - - // Then X1 EMF on the face is *averaged* along X3 - b = KDomain::GetRange(rc, domain, E1, coarse); - IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, E1, coarse); - const int jf = (binner) ? bi.js : bi.je; // j index of polar face - parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_EMF1_" + bname, pmb->exec_space, - 0, 1, b.is, b.ie, - KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int& i) { - // Sum the (non-ghost) X1 direction fluxes along the pole at zone i - // Recall both faces fall in the domain, so we neglect the right - double emf_sum; - Kokkos::Sum sum_reducer(emf_sum); - parthenon::par_reduce_inner(member, bi.ks, bi.ke - 1, - [&](const int& k, double& local_result) { - local_result += emfpack(E1, 0, k, jf, i); - } - , sum_reducer); - - // Calculate the average and set all EMFs identically (even ghosts, to keep divB) - const double emf_av = emf_sum / (bi.ke - bi.ks); - parthenon::par_for_inner(member, b.ks, b.ke, - [&](const int& k) { - emfpack(E1, 0, k, jf, i) = emf_av; - } - ); - } - ); -} - -void B_CT::DestructiveBoundaryClean(MeshBlockData *rc, IndexDomain domain, const VariablePack &fpack, bool coarse) -{ - // Set XN faces to keep clean divergence at outflow XN boundary (nearly always X1) - // Feels wrong to work backward from no divergence, but they are just outflow... - auto pmb = rc->GetBlockPointer(); - const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); - const std::string bname = KBoundaries::BoundaryName(bface); - const int bdir = KBoundaries::BoundaryDirection(bface); - const bool binner = KBoundaries::BoundaryIsInner(bface); - const TopologicalElement face = FaceOf(bdir); - // Correct last domain face, too - auto b = KDomain::GetRange(rc, domain, face, (binner) ? 0 : -1, (binner) ? 1 : 0, coarse); - // Need the coordinates for this boundary, uniquely - auto G = pmb->coords; - const int ndim = pmb->pmy_mesh->ndim; - if (domain == IndexDomain::inner_x1 || domain == IndexDomain::outer_x1) { - const int i_face = (binner) ? b.ie : b.is; - for (int iadd = 0; iadd <= (b.ie - b.is); iadd++) { - const int i = (binner) ? i_face - iadd : i_face + iadd; - const int last_rank_f = (binner) ? i + 1 : i - 1; - const int last_rank_c = (binner) ? i : i - 1; - const int outward_sign = (binner) ? -1. : 1.; - pmb->par_for( - "correct_face_vector_" + bname, b.ks, b.ke, b.js, b.je, i, i, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - // Other faces have been updated, just need to clean divergence - // Subtract off their contributions to find ours. Note our partner face contributes differently, - // depending on whether we're the i+1 "outward" face, or the i "innward" face - Real new_face = - (-outward_sign) * fpack(F1, 0, k, j, last_rank_f) * G.Volume(k, j, last_rank_f) - - (fpack(F2, 0, k, j + 1, last_rank_c) * G.Volume(k, j + 1, last_rank_c) - - fpack(F2, 0, k, j, last_rank_c) * G.Volume(k, j, last_rank_c)); - if (ndim > 2) - new_face -= fpack(F3, 0, k + 1, j, last_rank_c) * G.Volume(k + 1, j, last_rank_c) - - fpack(F3, 0, k, j, last_rank_c) * G.Volume(k, j, last_rank_c); - - fpack(F1, 0, k, j, i) = outward_sign * new_face / G.Volume(k, j, i); - } - ); - } - } else { - throw std::runtime_error("Divergence-free outflow replacement only implemented in X1!"); - } -} - double B_CT::MaxDivB(MeshData *md) { auto pmesh = md->GetMeshPointer(); @@ -591,7 +463,7 @@ double B_CT::MaxDivB(MeshData *md) pmb0->par_reduce("divB_max", block.s, block.e, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, double &local_result) { const auto& G = B_U.GetCoords(b); - double local_divb = face_div(G, B_U(b), ndim, k, j, i); + double local_divb = m::abs(face_div(G, B_U(b), ndim, k, j, i)); if (local_divb > local_result) local_result = local_divb; } , max_reducer); @@ -614,14 +486,13 @@ double B_CT::BlockMaxDivB(MeshBlockData *rc) pmb->par_reduce("divB_max", b.ks, b.ke, b.js, b.je, b.is, b.ie, KOKKOS_LAMBDA (const int &k, const int &j, const int &i, double &local_result) { const auto& G = B_U.GetCoords(); - double local_divb = face_div(G, B_U, ndim, k, j, i); + double local_divb = m::abs(face_div(G, B_U, ndim, k, j, i)); if (local_divb > local_result) local_result = local_divb; } , max_reducer); return max_divb; } - double B_CT::GlobalMaxDivB(MeshData *md, bool all_reduce) { if (all_reduce) { diff --git a/kharma/b_ct/b_ct.hpp b/kharma/b_ct/b_ct.hpp index 87defd8a..e9b1a605 100644 --- a/kharma/b_ct/b_ct.hpp +++ b/kharma/b_ct/b_ct.hpp @@ -122,12 +122,18 @@ void CalcDivB(MeshData *md, std::string divb_field_name="divB"); * * *mostly. I think */ -void ZeroEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse); +void ZeroBoundaryEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse); /** * Average all EMFs corresponding to the coordinate pole location, e.g. usually all E1 on X2 faces */ -void AverageEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse); +void AverageBoundaryEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse); + +/** + * Subtract the average B3 from each face, as if a loop reconnected across the polar boundary. + * Preserves divB, since differences across cells remain the same after subtracting a constant. + */ +void ReconnectBoundaryB3(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse); /** * Reset an outflow condition to have no divergence, even if a field line exits the domain. @@ -135,4 +141,22 @@ void AverageEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack< */ void DestructiveBoundaryClean(MeshBlockData *rc, IndexDomain domain, const VariablePack &fpack, bool coarse); +/** + * Take the curl over the whole domain. Defined in-header since it's templated on face and NDIM + */ +template +inline void EdgeCurl(MeshBlockData *rc, const GridVector& A, + const VariablePack& B_U, IndexDomain domain) +{ + auto pmb = rc->GetBlockPointer(); + const auto &G = pmb->coords; + IndexRange3 bB = KDomain::GetRange(rc, domain, el); + pmb->par_for( + "EdgeCurl", bB.ks, bB.ke, bB.js, bB.je, bB.is, bB.ie, + KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { + B_CT::edge_curl(G, A, B_U, k, j, i); + } + ); +} + } diff --git a/kharma/b_ct/b_ct_boundaries.cpp b/kharma/b_ct/b_ct_boundaries.cpp new file mode 100644 index 00000000..58f3639c --- /dev/null +++ b/kharma/b_ct/b_ct_boundaries.cpp @@ -0,0 +1,370 @@ +/* + * File: b_ct_boundaries.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 + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +#include "b_ct.hpp" + +#include "decs.hpp" +#include "domain.hpp" +#include "grmhd.hpp" +#include "grmhd_functions.hpp" +#include "kharma.hpp" + +void B_CT::ZeroBoundaryEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse) +{ + auto pmb = rc->GetBlockPointer(); + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const std::string bname = KBoundaries::BoundaryName(bface); + const int bdir = KBoundaries::BoundaryDirection(bface); + const bool binner = KBoundaries::BoundaryIsInner(bface); + // Select edges which lie on the domain face, zero only those + for (auto &el : OrthogonalEdges(bdir)) { + auto b = KDomain::GetBoundaryRange(rc, domain, el, coarse); + const int jf = (binner) ? b.je : b.js; + pmb->par_for( + "zero_EMF_" + bname, b.ks, b.ke, jf, jf, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + emfpack(el, 0, k, j, i) = 0; + } + ); + } +} + +void B_CT::AverageBoundaryEMF(MeshBlockData *rc, IndexDomain domain, const VariablePack &emfpack, bool coarse) +{ + auto pmb = rc->GetBlockPointer(); + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const std::string bname = KBoundaries::BoundaryName(bface); + const int bdir = KBoundaries::BoundaryDirection(bface); + const bool binner = KBoundaries::BoundaryIsInner(bface); + const int ndim = KDomain::GetNDim(rc); + + for (auto &el : OrthogonalEdges(bdir)) { + if (bdir == X2DIR && el == E3 && pmb->coords.coords.is_spherical()) { + // X3 EMF must be zero *on* polar face, since edge size is 0 + IndexRange3 b = KDomain::GetBoundaryRange(rc, domain, el, coarse); + pmb->par_for( + "zero_polar_EMF3_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + emfpack(el, 0, k, j, i) = 0; + } + ); + } else if (ndim < 3 && ((bdir == X2DIR && el == E1) || (bdir == X1DIR && el == E2))) { + // In 2D, "averaging" should just mean not zeroing E1 on X2 or E2 on X1 + continue; + } else { + // Otherwise the EMF at `el` is *averaged* along its perpendicular direction + IndexRange3 b = KDomain::GetRange(rc, domain, el, coarse); + IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, el, coarse); + // Calculate face index and outer sum index + int cface, inner_dir; + IndexRange outer; + if (bdir == X1DIR) { + cface = (binner) ? bi.is : bi.ie; + if (el == E2) { + outer = {b.js, b.je}; + inner_dir = X3DIR; + } else { + outer = {b.ks, b.ke}; + inner_dir = X2DIR; + } + } else if (bdir == X2DIR) { + cface = (binner) ? bi.js : bi.je; + if (el == E1) { // COMMON CASE + outer = {b.is, b.ie}; + inner_dir = X3DIR; + } else { + outer = {b.ks, b.ke}; + inner_dir = X1DIR; + } + } else { + cface = (binner) ? bi.ks : bi.ke; + if (el == E1) { + outer = {b.is, b.ie}; + inner_dir = X2DIR; + } else { + outer = {b.js, b.je}; + inner_dir = X1DIR; + } + } + parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_EMF1_" + bname, pmb->exec_space, + 0, 1, outer.s, outer.e, + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int& o) { + double emf_sum = 0.; + Kokkos::Sum sum_reducer(emf_sum); + + // One of these won't be used + // Outer loop is along face corresponding to our element + const int ii = (el == E1) ? o : cface; + const int jj = (el == E2) ? o : cface; + const int kk = (el == E3) ? o : cface; + + // Sum the non-ghost fluxes in our desired averaging direction + int len; + if (inner_dir == X1DIR) { + len = bi.ie - bi.is; + parthenon::par_reduce_inner(member, bi.is, bi.ie - 1, + [&](const int& i, double& local_result) { + local_result += emfpack(el, 0, kk, jj, i); + } + , sum_reducer); + } else if (inner_dir == X2DIR) { + len = bi.je - bi.js; + parthenon::par_reduce_inner(member, bi.js, bi.je - 1, + [&](const int& j, double& local_result) { + local_result += emfpack(el, 0, kk, j, ii); + } + , sum_reducer); + } else { + len = bi.ke - bi.ks; + parthenon::par_reduce_inner(member, bi.ks, bi.ke - 1, + [&](const int& k, double& local_result) { + local_result += emfpack(el, 0, k, jj, ii); + } + , sum_reducer); + } + + // Calculate the average + const double emf_av = emf_sum / len; + + // Set all EMFs identically (even ghosts, to keep divB) + if (inner_dir == X1DIR) { + parthenon::par_for_inner(member, b.is, b.ie, + [&](const int& i) { + emfpack(el, 0, kk, jj, i) = emf_av; + } + ); + } else if (inner_dir == X2DIR) { + parthenon::par_for_inner(member, b.js, b.je, + [&](const int& j) { + emfpack(el, 0, kk, j, ii) = emf_av; + } + ); + } else { + parthenon::par_for_inner(member, b.ks, b.ke, + [&](const int& k) { + emfpack(el, 0, k, jj, ii) = emf_av; + } + ); + } + } + ); + } + } +} + +void B_CT::DestructiveBoundaryClean(MeshBlockData *rc, IndexDomain domain, const VariablePack &fpack, bool coarse) +{ + // Set XN faces to keep clean divergence at outflow XN boundary + // Feels wrong to work backward from no divergence, but they are just outflow... + auto pmb = rc->GetBlockPointer(); + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const std::string bname = KBoundaries::BoundaryName(bface); + const int bdir = KBoundaries::BoundaryDirection(bface); + const bool binner = KBoundaries::BoundaryIsInner(bface); + const TopologicalElement face = FaceOf(bdir); + // Correct last domain face, too + auto b = KDomain::GetRange(rc, domain, face, (binner) ? 0 : -1, (binner) ? 1 : 0, coarse); + // Need the coordinates for this boundary, uniquely + auto G = pmb->coords; + const int ndim = pmb->pmy_mesh->ndim; + if (bdir == X1DIR) { + const int i_face = (binner) ? b.ie : b.is; + for (int iadd = 0; iadd <= (b.ie - b.is); iadd++) { + const int i = (binner) ? i_face - iadd : i_face + iadd; + const int last_rank_f = (binner) ? i + 1 : i - 1; + const int last_rank_c = (binner) ? i : i - 1; + const int outward_sign = (binner) ? -1. : 1.; + pmb->par_for( + "correct_face_vector_" + bname, b.ks, b.ke, b.js, b.je, i, i, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + // Other faces have been updated, just need to clean divergence + // Subtract off their contributions to find ours. Note our partner face contributes differently, + // depending on whether we're the i+1 "outward" face, or the i "innward" face + Real new_face = - (-outward_sign) * fpack(F1, 0, k, j, last_rank_f) * G.Volume(k, j, last_rank_f) + - (fpack(F2, 0, k, j + 1, last_rank_c) * G.Volume(k, j + 1, last_rank_c) + - fpack(F2, 0, k, j, last_rank_c) * G.Volume(k, j, last_rank_c)); + if (ndim > 2) + new_face -= fpack(F3, 0, k + 1, j, last_rank_c) * G.Volume(k + 1, j, last_rank_c) + - fpack(F3, 0, k, j, last_rank_c) * G.Volume(k, j, last_rank_c); + + fpack(F1, 0, k, j, i) = outward_sign * new_face / G.Volume(k, j, i); + } + ); + } + } else if (bdir == X2DIR) { + const int j_face = (binner) ? b.je : b.js; + for (int jadd = 0; jadd <= (b.je - b.js); jadd++) { + const int j = (binner) ? j_face - jadd : j_face + jadd; + const int last_rank_f = (binner) ? j + 1 : j - 1; + const int last_rank_c = (binner) ? j : j - 1; + const int outward_sign = (binner) ? -1. : 1.; + pmb->par_for( + "correct_face_vector_" + bname, b.ks, b.ke, j, j, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + Real new_face = - (-outward_sign) * fpack(F2, 0, k, last_rank_f, i) * G.Volume(k, last_rank_f, i) + - (fpack(F1, 0, k, last_rank_c, i + 1) * G.Volume(k, last_rank_c, i + 1) + - fpack(F1, 0, k, last_rank_c, i) * G.Volume(k, last_rank_c, i)); + if (ndim > 2) + new_face -= fpack(F3, 0, k + 1, last_rank_c, i) * G.Volume(k + 1, last_rank_c, i) + - fpack(F3, 0, k, last_rank_c, i) * G.Volume(k, j, last_rank_c, i); + + fpack(F2, 0, k, j, i) = outward_sign * new_face / G.Volume(k, j, i); + } + ); + } + } else { + const int k_face = (binner) ? b.ie : b.is; + for (int kadd = 0; kadd <= (b.ie - b.is); kadd++) { + const int k = (binner) ? k_face - kadd : k_face + kadd; + const int last_rank_f = (binner) ? k + 1 : k - 1; + const int last_rank_c = (binner) ? k : k - 1; + const int outward_sign = (binner) ? -1. : 1.; + pmb->par_for( + "correct_face_vector_" + bname, k, k, b.js, b.je, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + Real new_face = - (-outward_sign) * fpack(F3, 0, last_rank_f, j, i) * G.Volume(last_rank_f, j, i) + - (fpack(F1, 0, last_rank_c, j, i + 1) * G.Volume(last_rank_c, j, i + 1) + - fpack(F1, 0, last_rank_c, j, i) * G.Volume(last_rank_c, j, i)); + - (fpack(F2, 0, last_rank_c, j + 1, i) * G.Volume(last_rank_c, j + 1, i) + - fpack(F2, 0, last_rank_c, j, i) * G.Volume(last_rank_c, j, i)); + + fpack(F3, 0, k, j, i) = outward_sign * new_face / G.Volume(k, j, i); + } + ); + } + } +} + +// Reducer class for unused MinAbs B below +// template +// struct MinAbsReducer { +// public: +// // Required +// typedef MinAbsReducer reducer; +// typedef double value_type; +// // typedef Kokkos::View +// // result_view_type; + +// private: +// value_type& value; + +// public: +// KOKKOS_INLINE_FUNCTION +// MinAbsReducer(value_type& value_) : value(value_) {} + +// // Required +// KOKKOS_INLINE_FUNCTION +// void join(value_type& dest, const value_type& src) const { +// dest = (m::abs(src) < m::abs(dest)) ? src : dest; +// } + +// KOKKOS_INLINE_FUNCTION +// void init(value_type& val) const { val = 0; } + +// KOKKOS_INLINE_FUNCTION +// value_type& reference() const { return value; } + +// //KOKKOS_INLINE_FUNCTION +// //result_view_type view() const { return result_view_type(&value, 1); } + +// KOKKOS_INLINE_FUNCTION +// bool references_scalar() const { return true; } +// }; + +void B_CT::ReconnectBoundaryB3(MeshBlockData *rc, IndexDomain domain, const VariablePack &fpack, bool coarse) +{ + // We're also sometimes called on coarse buffers with or without AMR. + // Use of transmitting polar conditions when coarse buffers matter (e.g., refinement + // boundary touching the pole) is UNSUPPORTED + if (coarse) return; + + // Pull boundary properties + auto pmb = rc->GetBlockPointer(); + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const bool binner = KBoundaries::BoundaryIsInner(bface); + const int bdir = KBoundaries::BoundaryDirection(bface); + const auto bname = KBoundaries::BoundaryName(bface); + + // Subtract the average B3 as "reconnection" + IndexRange3 b = KDomain::GetRange(rc, domain, F3, coarse); + IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, F3, coarse); + const int jf = (binner) ? bi.js : bi.je; // j index of last zone next to pole + parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_B3_" + bname, pmb->exec_space, + 0, 1, 0, fpack.GetDim(4)-1, b.is, b.ie, + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int &v, const int& i) { + // Sum the first rank of B3 + double B3_sum = 0.; + Kokkos::Sum sum_reducer(B3_sum); + parthenon::par_reduce_inner(member, bi.ks, bi.ke - 1, + [&](const int& k, double& local_result) { + local_result += fpack(F3, v, k, jf, i); + } + , sum_reducer); + + // Calculate the average and modify all B3 identically + // This will preserve their differences->divergence + const double B3_av = B3_sum / (bi.ke - bi.ks); + parthenon::par_for_inner(member, b.ks, b.ke, + [&](const int& k) { + fpack(F3, v, k, jf, i) -= B3_av; + } + ); + } + ); + // Option for subtracting minimum by absolute value, much less stable + // parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_B3_" + bname, pmb->exec_space, + // 0, 1, 0, fpack.GetDim(4)-1, b.is, b.ie, + // KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int &v, const int& i) { + // // Sum the first rank of B3 + // double B3_min = 0.; + // MinAbsReducer min_abs_reducer(B3_min); + // parthenon::par_reduce_inner(member, bi.ks, bi.ke - 1, + // [&](const int& k, double& local_result) { + // // Compare unsigned + // if (m::abs(fpack(F3, v, k, jf, i)) < m::abs(local_result)) { + // // Assign signed, reducer will compare unsigned + // local_result = fpack(F3, v, k, jf, i); + // } + // } + // , min_abs_reducer); + + // // Subtract from all B3 identically + // // This will preserve their differences->divergence + // parthenon::par_for_inner(member, b.ks, b.ke, + // [&](const int& k) { + // fpack(F3, v, k, jf, i) -= B3_min; + // } + // ); + // } + // ); +} \ No newline at end of file diff --git a/kharma/b_ct/b_ct_functions.hpp b/kharma/b_ct/b_ct_functions.hpp index d8d9072e..a294d0fe 100644 --- a/kharma/b_ct/b_ct_functions.hpp +++ b/kharma/b_ct/b_ct_functions.hpp @@ -100,28 +100,13 @@ KOKKOS_INLINE_FUNCTION void edge_curl(const GRCoordinates& G, const GridVector& } } -template -KOKKOS_INLINE_FUNCTION void EdgeCurl(MeshBlockData *rc, const GridVector& A, - const VariablePack& B_U, IndexDomain domain) -{ - auto pmb = rc->GetBlockPointer(); - const auto &G = pmb->coords; - IndexRange3 bB = KDomain::GetRange(rc, domain, el); - pmb->par_for( - "EdgeCurl", bB.ks, bB.ke, bB.js, bB.je, bB.is, bB.ie, - KOKKOS_LAMBDA(const int &k, const int &j, const int &i) { - B_CT::edge_curl(G, A, B_U, k, j, i); - } - ); -} - KOKKOS_INLINE_FUNCTION Real upwind_diff(const VariableFluxPack& B_U, const VariablePack& emfc, const VariablePack& uvec, const int& comp, const int& dir, const int& vdir, const int& k, const int& j, const int& i, const bool& left_deriv) { // See SG09 eq 23 // Upwind based on vel(vdir) at the left face in vdir (contact mode) - TopologicalElement face = FaceOf(vdir); + TopologicalElement face = (vdir == 1) ? F1 : ((vdir == 2) ? F2 : F3); // const Real contact_vel = uvec(face, vdir-1, k, j, i); // Upwind by one zone in dir const int i_up = (vdir == 1) ? i - 1 : i; @@ -157,6 +142,7 @@ KOKKOS_INLINE_FUNCTION Real upwind_diff(const VariableFluxPack& B_U, const } // Only through formatting has the following been made even a little comprehensible. +// TODO move it into Parthenon! template KOKKOS_FORCEINLINE_FUNCTION Real F(const ParArrayND &fine, const Coordinates_t &coords, int l, int m, int n, int fk, int fj, int fi) diff --git a/kharma/b_flux_ct/b_flux_ct.cpp b/kharma/b_flux_ct/b_flux_ct.cpp index 23176810..62a1d717 100644 --- a/kharma/b_flux_ct/b_flux_ct.cpp +++ b/kharma/b_flux_ct/b_flux_ct.cpp @@ -37,6 +37,7 @@ #include "b_flux_ct.hpp" #include "decs.hpp" +#include "domain.hpp" #include "grmhd.hpp" #include "kharma.hpp" @@ -53,7 +54,7 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

GetBoolean("coordinates", "spherical"); bool fix_polar_flux = pin->GetOrAddBoolean("b_field", "fix_polar_flux", spherical); params.Add("fix_polar_flux", fix_polar_flux); @@ -61,17 +62,23 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

GetOrAddBoolean("b_field", "fix_flux_x1", false); - // Split out options. Turn off inner edge by default if inner bound is within EH + // Split out options. Turn off inner edge by default if inner bound is within EH (TODO turn on/off based on BOUNDARIES) bool r_in_eh = spherical && pin->GetBoolean("coordinates", "domain_intersects_eh"); bool fix_flux_inner_x1 = pin->GetOrAddBoolean("b_field", "fix_flux_inner_x1", fix_flux_x1 && !r_in_eh); params.Add("fix_flux_inner_x1", fix_flux_inner_x1); bool fix_flux_outer_x1 = pin->GetOrAddBoolean("b_field", "fix_flux_outer_x1", fix_flux_x1); params.Add("fix_flux_outer_x1", fix_flux_outer_x1); + // Bflux0 on x2. NOT good for polar boundaries + bool fix_flux_x2 = pin->GetOrAddBoolean("b_field", "fix_flux_x2", false); + bool fix_flux_inner_x2 = pin->GetOrAddBoolean("b_field", "fix_flux_inner_x2", fix_flux_x2); + params.Add("fix_flux_inner_x2", fix_flux_inner_x2); + bool fix_flux_outer_x2 = pin->GetOrAddBoolean("b_field", "fix_flux_outer_x2", fix_flux_x2); + params.Add("fix_flux_outer_x2", fix_flux_outer_x2); // This reverts to a more ham-fisted fix which explicitly disallows flux crossing the X1 face. // This version requires *inverted* B1 across the face, potentially just using reflecting conditions for B // Using this version is tremendously inadvisable: consult your simulator before applying. - bool use_old_x1_fix = pin->GetOrAddBoolean("b_field", "use_old_x1_fix", false); - params.Add("use_old_x1_fix", use_old_x1_fix); + bool use_old_flux_fix = pin->GetOrAddBoolean("b_field", "use_old_flux_fix", false); + params.Add("use_old_flux_fix", use_old_flux_fix); // KHARMA requires some kind of field transport if there is a magnetic field allocated // Use this if you actually want to disable all magnetic field flux corrections, @@ -255,18 +262,41 @@ void BlockPtoU(MeshBlockData *rc, IndexDomain domain, bool coarse) void FixFlux(MeshData *md) { - // TODO flags here auto pmb0 = md->GetBlockData(0)->GetBlockPointer(); auto& params = pmb0->packages.Get("B_FluxCT")->AllParams(); + // Poles specifically can't use Bflux0 if (params.Get("fix_polar_flux")) { - FixBoundaryFlux(md, IndexDomain::inner_x2, false); - FixBoundaryFlux(md, IndexDomain::outer_x2, false); + ZeroBoundaryFlux(md, IndexDomain::inner_x2, false); + ZeroBoundaryFlux(md, IndexDomain::outer_x2, false); } + // Everything else should default to Bflux0 if (params.Get("fix_flux_inner_x1")) { - FixBoundaryFlux(md, IndexDomain::inner_x1, false); + if (params.Get("use_old_flux_fix")) { + ZeroBoundaryFlux(md, IndexDomain::inner_x1, false); + } else { + Bflux0(md, IndexDomain::inner_x1, false); + } } if (params.Get("fix_flux_outer_x1")) { - FixBoundaryFlux(md, IndexDomain::outer_x1, false); + if (params.Get("use_old_flux_fix")) { + ZeroBoundaryFlux(md, IndexDomain::outer_x1, false); + } else { + Bflux0(md, IndexDomain::outer_x1, false); + } + } + if (params.Get("fix_flux_inner_x2")) { + if (params.Get("use_old_flux_fix")) { + ZeroBoundaryFlux(md, IndexDomain::inner_x2, false); + } else { + Bflux0(md, IndexDomain::inner_x2, false); + } + } + if (params.Get("fix_flux_outer_x2")) { + if (params.Get("use_old_flux_fix")) { + ZeroBoundaryFlux(md, IndexDomain::inner_x2, false); + } else { + Bflux0(md, IndexDomain::outer_x2, false); + } } FluxCT(md); } @@ -338,21 +368,14 @@ void FluxCT(MeshData *md) } } -void FixBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) +void ZeroBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) { + // TODO write ONE implementation for any boundary auto pmesh = md->GetMeshPointer(); auto pmb0 = pmesh->block_list[0]; const int ndim = pmesh->ndim; if (ndim < 2) return; - // Option for old, pre-Bflux0 - const bool use_old_x1_fix = pmb0->packages.Get("B_FluxCT")->Param("use_old_x1_fix"); - - auto bounds = coarse ? pmb0->c_cellbounds : pmb0->cellbounds; - const IndexRange ib = bounds.GetBoundsI(IndexDomain::interior); - const IndexRange jb = bounds.GetBoundsJ(IndexDomain::interior); - const IndexRange kb = bounds.GetBoundsK(IndexDomain::interior); - // Imagine a corner of the domain, with ghost and physical zones // as below, denoted w/'g' and 'p' respectively. // ... @@ -368,7 +391,10 @@ void FixBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) // Therefore in e.g. X1 faces, we need to update fluxes on the domain: // [0,N1+1],[-1,N2+1],[-1,N3+1] // These indices arrange for that. - + auto bounds = coarse ? pmb0->c_cellbounds : pmb0->cellbounds; + const IndexRange ib = bounds.GetBoundsI(IndexDomain::interior); + const IndexRange jb = bounds.GetBoundsJ(IndexDomain::interior); + const IndexRange kb = bounds.GetBoundsK(IndexDomain::interior); // For faces const IndexRange ibf = IndexRange{ib.s, ib.e + 1}; const IndexRange jbf = IndexRange{jb.s, jb.e + 1}; @@ -392,8 +418,8 @@ void FixBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { B_F.flux(X2DIR, V1, k, j, i) = 0.; B_F.flux(X2DIR, V3, k, j, i) = 0.; - B_F.flux(X1DIR, V2, k, j - 1, i) = -B_F.flux(X1DIR, V2, k, j, i); - if (ndim > 2) B_F.flux(X3DIR, V2, k, j - 1, i) = -B_F.flux(X3DIR, V2, k, j, i); + B_F.flux(X1DIR, V2, k, j-1, i) = -B_F.flux(X1DIR, V2, k, j, i); + if (ndim > 2) B_F.flux(X3DIR, V2, k, j-1, i) = -B_F.flux(X3DIR, V2, k, j, i); } ); } @@ -409,70 +435,118 @@ void FixBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse) } ); } + // These boundary conditions need to arrange for B1 to be inverted in ghost cells. + // This is no longer pure outflow, but might be thought of as a "nicer" version of + // reflecting conditions: + // 1. Since B1 is inverted, B1 on the domain face will tend to 0 (it's not quite reflected, but basically) + // (obviously don't enable this for monopole test problems!) + // 2. However, B2 and B3 are normal outflow conditions -- despite the fluxes here, the outflow + // conditions will set them equal to the last zone. + if (domain == IndexDomain::inner_x1 && + pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) { + pmb->par_for("fix_flux_b_in_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + B_F.flux(X1DIR, V2, k, j, i) = 0.; + B_F.flux(X1DIR, V3, k, j, i) = 0.; + B_F.flux(X2DIR, V1, k, j, i - 1) = -B_F.flux(X2DIR, V1, k, j, i); + if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i - 1) = -B_F.flux(X3DIR, V1, k, j, i); + } + ); + } - // TODO(BSP) could check here we're operating with the right boundaries: Dirichlet for Bflux0, - // reflecting/B1 reflect for old stuff - if (!use_old_x1_fix) { - // "Bflux0" prescription for keeping divB~=0 on zone corners of the interior & exterior X1 faces - // Courtesy of & implemented by Hyerin Cho - // Allows nonzero flux across X1 boundary but still keeps divB=0 (turns out effectively to have 0 flux) - // Usable only for Dirichlet conditions - if (domain == IndexDomain::inner_x1 && - pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) - { - pmb->par_for("fix_flux_b_in", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s, // Hyerin (12/28/22) for 1st & 2nd prescription - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - // Allows nonzero flux across X1 boundary but still keeps divB=0 (turns out effectively to have 0 flux) - if (ndim > 1) B_F.flux(X2DIR, V1, k, j, i-1) = -B_F.flux(X2DIR, V1, k, j, i) + B_F.flux(X1DIR, V2, k, j, i) + B_F.flux(X1DIR, V2, k, j-1, i); - if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i-1) = -B_F.flux(X3DIR, V1, k, j, i) + B_F.flux(X1DIR, V3, k, j, i) + B_F.flux(X1DIR, V3, k-1, j, i); - } - ); + if (domain == IndexDomain::outer_x1 && + pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) { + pmb->par_for("fix_flux_b_out_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + B_F.flux(X1DIR, V2, k, j, i) = 0.; + B_F.flux(X1DIR, V3, k, j, i) = 0.; + B_F.flux(X2DIR, V1, k, j, i) = -B_F.flux(X2DIR, V1, k, j, i - 1); + if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i) = -B_F.flux(X3DIR, V1, k, j, i - 1); + } + ); + } + } +} - } - if (domain == IndexDomain::outer_x1 && - pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) - { - pmb->par_for("fix_flux_b_out", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e, // Hyerin (12/28/22) for 1st & 2nd prescription - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - // (02/06/23) 2nd prescription that allows nonzero flux across X1 boundary but still keeps divB=0 - if (ndim > 1) B_F.flux(X2DIR, V1, k, j, i) = -B_F.flux(X2DIR, V1, k, j, i-1) + B_F.flux(X1DIR, V2, k, j, i) + B_F.flux(X1DIR, V2, k, j-1, i); - if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i) = -B_F.flux(X3DIR, V1, k, j, i-1) + B_F.flux(X1DIR, V3, k, j, i) + B_F.flux(X1DIR, V3, k-1, j, i); - } - ); - } - } else { - // These boundary conditions need to arrange for B1 to be inverted in ghost cells. - // This is no longer pure outflow, but might be thought of as a "nicer" version of - // reflecting conditions: - // 1. Since B1 is inverted, B1 on the domain face will tend to 0 (it's not quite reflected, but basically) - // (obviously don't enable this for monopole test problems!) - // 2. However, B2 and B3 are normal outflow conditions -- despite the fluxes here, the outflow - // conditions will set them equal to the last zone. - if (domain == IndexDomain::inner_x1 && - pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) { - pmb->par_for("fix_flux_b_in_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - B_F.flux(X1DIR, V2, k, j, i) = 0.; - B_F.flux(X1DIR, V3, k, j, i) = 0.; - B_F.flux(X2DIR, V1, k, j, i - 1) = -B_F.flux(X2DIR, V1, k, j, i); - if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i - 1) = -B_F.flux(X3DIR, V1, k, j, i); - } - ); - } +void Bflux0(MeshData *md, IndexDomain domain, bool coarse) +{ + auto pmesh = md->GetMeshPointer(); + auto pmb0 = pmesh->block_list[0]; + const int ndim = pmesh->ndim; + if (ndim < 2) return; + + // Indices, see ZeroBoundaryFlux + auto bounds = coarse ? pmb0->c_cellbounds : pmb0->cellbounds; + const IndexRange ib = bounds.GetBoundsI(IndexDomain::interior); + const IndexRange jb = bounds.GetBoundsJ(IndexDomain::interior); + const IndexRange kb = bounds.GetBoundsK(IndexDomain::interior); + // For faces + const IndexRange ibf = IndexRange{ib.s, ib.e + 1}; + const IndexRange jbf = IndexRange{jb.s, jb.e + 1}; + // Won't need X3 faces + //const IndexRange kbf = IndexRange{kb.s, kb.e + (ndim > 2)}; + // For sides + const IndexRange ibs = IndexRange{ib.s - 1, ib.e + 1}; + const IndexRange jbs = IndexRange{jb.s - (ndim > 1), jb.e + (ndim > 1)}; + const IndexRange kbs = IndexRange{kb.s - (ndim > 2), kb.e + (ndim > 2)}; + + // Make sure the polar EMFs are 0 when performing fluxCT + // 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& B_F = rc->PackVariablesAndFluxes(std::vector{"cons.B"}); + + // "Bflux0" prescription for keeping divB~=0 on zone corners of the interior & exterior X1 faces + // Courtesy of & implemented by Hyerin Cho + // Allows nonzero flux across X1 boundary but still keeps divB=0 (turns out effectively to have 0 flux) + // Usable only for Dirichlet conditions + if (domain == IndexDomain::inner_x1 && + pmb->boundary_flag[BoundaryFace::inner_x1] == BoundaryFlag::user) + { + pmb->par_for("fix_flux_b_in", kbs.s, kbs.e, jbs.s, jbs.e, ibf.s, ibf.s, // Hyerin (12/28/22) for 1st & 2nd prescription + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + // Allows nonzero flux across X1 boundary but still keeps divB=0 (turns out effectively to have 0 flux) + if (ndim > 1) B_F.flux(X2DIR, V1, k, j, i-1) = -B_F.flux(X2DIR, V1, k, j, i) + B_F.flux(X1DIR, V2, k, j, i) + B_F.flux(X1DIR, V2, k, j-1, i); + if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i-1) = -B_F.flux(X3DIR, V1, k, j, i) + B_F.flux(X1DIR, V3, k, j, i) + B_F.flux(X1DIR, V3, k-1, j, i); + } + ); + + } + if (domain == IndexDomain::inner_x2 && + pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) + { + pmb->par_for("fix_flux_b_in", kbs.s, kbs.e, jbf.s, jbf.s, ibs.s, ibs.e, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + if (ndim > 2) B_F.flux(X3DIR, V2, k, j-1, i) = -B_F.flux(X3DIR, V2, k, j, i) + B_F.flux(X2DIR, V3, k, j, i) + B_F.flux(X2DIR, V3, k-1, j, i); + if (ndim > 1) B_F.flux(X1DIR, V2, k, j-1, i) = -B_F.flux(X1DIR, V2, k, j, i) + B_F.flux(X2DIR, V1, k, j, i) + B_F.flux(X2DIR, V1, k, j, i-1); + } + ); - if (domain == IndexDomain::outer_x1 && - pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) { - pmb->par_for("fix_flux_b_out_old", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - B_F.flux(X1DIR, V2, k, j, i) = 0.; - B_F.flux(X1DIR, V3, k, j, i) = 0.; - B_F.flux(X2DIR, V1, k, j, i) = -B_F.flux(X2DIR, V1, k, j, i - 1); - if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i) = -B_F.flux(X3DIR, V1, k, j, i - 1); - } - ); - } } + // OUTER + if (domain == IndexDomain::outer_x1 && + pmb->boundary_flag[BoundaryFace::outer_x1] == BoundaryFlag::user) + { + pmb->par_for("fix_flux_b_out", kbs.s, kbs.e, jbs.s, jbs.e, ibf.e, ibf.e, // Hyerin (12/28/22) for 1st & 2nd prescription + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + // (02/06/23) 2nd prescription that allows nonzero flux across X1 boundary but still keeps divB=0 + if (ndim > 1) B_F.flux(X2DIR, V1, k, j, i) = -B_F.flux(X2DIR, V1, k, j, i-1) + B_F.flux(X1DIR, V2, k, j, i) + B_F.flux(X1DIR, V2, k, j-1, i); + if (ndim > 2) B_F.flux(X3DIR, V1, k, j, i) = -B_F.flux(X3DIR, V1, k, j, i-1) + B_F.flux(X1DIR, V3, k, j, i) + B_F.flux(X1DIR, V3, k-1, j, i); + } + ); + } + if (domain == IndexDomain::outer_x2 && + pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) + { + pmb->par_for("fix_flux_b_out", kbs.s, kbs.e, jbf.e, jbf.e, ibs.s, ibs.e, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + if (ndim > 2) B_F.flux(X3DIR, V2, k, j, i) = -B_F.flux(X3DIR, V2, k, j-1, i) + B_F.flux(X2DIR, V3, k, j, i) + B_F.flux(X2DIR, V3, k-1, j, i); + if (ndim > 1) B_F.flux(X1DIR, V2, k, j, i) = -B_F.flux(X1DIR, V2, k, j-1, i) + B_F.flux(X2DIR, V1, k, j, i) + B_F.flux(X2DIR, V1, k, j, i-1); + } + ); + } } } diff --git a/kharma/b_flux_ct/b_flux_ct.hpp b/kharma/b_flux_ct/b_flux_ct.hpp index 9c1700a0..8ea469a6 100644 --- a/kharma/b_flux_ct/b_flux_ct.hpp +++ b/kharma/b_flux_ct/b_flux_ct.hpp @@ -44,13 +44,13 @@ * * This physics package implements B field transport with Flux-CT (Toth 2000) * - * This requires only the values at cell centers + * This requires only the magnetic field value at cell centers * * This implementation includes conversion from "primitive" to "conserved" B and back */ namespace B_FluxCT { /** - * Declare fields, initialize (few) parameters + * Declare fields, initialize parameters */ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr& packages); @@ -74,18 +74,26 @@ void BlockPtoU(MeshBlockData *md, IndexDomain domain, bool coarse=false); void MeshPtoU(MeshData *md, IndexDomain domain, bool coarse=false); /** - * All flux corrections required by this package + * Apply all flux corrections required by this package to ensure small divB */ void FixFlux(MeshData *md); + /** * Modify the B field fluxes to take a constrained-transport step as in Toth (2000) */ void FluxCT(MeshData *md); + /** * Modify the B field fluxes just beyond the polar (or radial) boundary so as to * ensure no flux through the boundary after applying FluxCT */ -void FixBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse); +void ZeroBoundaryFlux(MeshData *md, IndexDomain domain, bool coarse); + +/** + * Modify the B field fluxes just beyond the radial (or polar) boundary so as to + * ensure the magnetic divergence is zero, even + */ +void Bflux0(MeshData *md, IndexDomain domain, bool coarse); /** * Alternate B field fix for X1 boundary, keeps zero divergence while permitting flux diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp index ba1741a5..aeb04782 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -33,6 +33,7 @@ */ #include "boundaries.hpp" +#include "bondi.hpp" #include "decs.hpp" #include "domain.hpp" #include "kharma.hpp" @@ -168,22 +169,31 @@ std::shared_ptr KBoundaries::Initialize(ParameterInput *pin, std: bool outflow_EMHD = pin->GetOrAddBoolean("boundaries", "outflow_EMHD_" + bname, false); params.Add("outflow_EMHD_" + bname, outflow_EMHD); - // Invert X2 face values to reflect across polar boundary - bool invert_F2 = pin->GetOrAddBoolean("boundaries", "reflect_face_vector_" + bname, (btype == "reflecting")); - params.Add("reflect_face_vector_"+bname, invert_F2); - // If you'll have field loops exiting the domain, outflow conditions need to be cleaned so as not to - // 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); - - // Special EMF averaging. Probably slow but beneficial for transmitting boundaries - bool average_EMF = pin->GetOrAddBoolean("boundaries", "average_EMF_" + bname, (btype == "transmitting")); - params.Add("average_EMF_"+bname, average_EMF); - // Otherwise, always zero EMFs to prevent B field escaping the domain in polar/dirichlet bounds - bool zero_EMF = pin->GetOrAddBoolean("boundaries", "zero_EMF_" + bname, ((bdir == X2DIR && spherical) - || (btype == "dirichlet")) - && !average_EMF); - params.Add("zero_EMF_"+bname, zero_EMF); + // Options specific to face-centered B fields, which require a lot of care at boundaries + if (packages->AllPackages().count("B_CT")) { + // Invert X2 face values to reflect across polar boundary + bool invert_F2 = pin->GetOrAddBoolean("boundaries", "reflect_face_vector_" + bname, (btype == "reflecting")); + params.Add("reflect_face_vector_"+bname, invert_F2); + // If you'll have field loops exiting the domain, outflow conditions need to be cleaned so as not to + // 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 pole w/face-CT. Maybe useful for reflecting too? + bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, (btype == "transmitting")); + params.Add("reconnect_B3_"+bname, reconnect_B3); + + // Special EMF averaging. Allows B slippage, e.g. around pole for transmitting conditions + // Still problems with averaging+dirichlet, maybe corners? + bool average_EMF = pin->GetOrAddBoolean("boundaries", "average_EMF_" + bname, (btype == "transmitting")); + params.Add("average_EMF_"+bname, average_EMF); + // Otherwise, always zero EMFs to prevent B field escaping the domain in polar/dirichlet bounds + bool zero_EMF = pin->GetOrAddBoolean("boundaries", "zero_EMF_" + bname, !average_EMF); + params.Add("zero_EMF_"+bname, zero_EMF); + } + // Advect together/cancel U3, under the theory it's in a similar position to B3 above (albeit no CT constraining it) + // Not enabled by default as it does not conserve angular momentum and isn't necessary for stability + bool cancel_U3 = pin->GetOrAddBoolean("boundaries", "cancel_U3_" + bname, false); + params.Add("cancel_U3_"+bname, cancel_U3); // String manip to get the Parthenon boundary name, e.g., "ox1_bc" @@ -300,6 +310,31 @@ std::shared_ptr KBoundaries::Initialize(ParameterInput *pin, std: default: break; } + } else if (btype == "bondi") { + // Boundaries will need these to be recorded into a 'Params' + AddBondiParameters(pin, *packages); + switch (bface) { + case BoundaryFace::inner_x1: + pkg->KBoundaries[bface] = SetBondi; + break; + case BoundaryFace::outer_x1: + pkg->KBoundaries[bface] = SetBondi; + break; + case BoundaryFace::inner_x2: + pkg->KBoundaries[bface] = SetBondi; + break; + case BoundaryFace::outer_x2: + pkg->KBoundaries[bface] = SetBondi; + break; + case BoundaryFace::inner_x3: + pkg->KBoundaries[bface] = SetBondi; + break; + case BoundaryFace::outer_x3: + pkg->KBoundaries[bface] = SetBondi; + break; + default: + break; + } } else { throw std::runtime_error("Unknown boundary type: "+btype); } @@ -337,6 +372,30 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD const auto bdir = BoundaryDirection(bface); const bool binner = BoundaryIsInner(bface); + // We get called over a lot of different packs depending on doing physical boundaries, EMF boundaries, + // boundaries during solves/GMG, and so on. Check once whether this is a "normal" boundary condition + // with the GRMHD variables + // TODO redesign boundary functions as per-package callbacks which take a pack+map + // TODO probably retire PackMHDPrims, it's not more useful than just packing on the flag. + PackIndexMap dummy_map; + bool full_grmhd_boundary = GRMHD::PackMHDPrims(rc.get(), dummy_map).GetDim(4) > 0; + + // Averaging ops on *physical* cells must be done before computing boundaries + // We should do a PreBoundaries callback... + if (pmb->packages.AllPackages().count("B_CT")) { + auto bfpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost, Metadata::GetUserFlag("B_CT")}); + if (params.Get("reconnect_B3_" + bname) && bfpack.GetDim(4) > 0) { + Flag("ReconnectFaceB_"+bname); + B_CT::ReconnectBoundaryB3(rc.get(), domain, bfpack, coarse); + EndFlag(); + } + } + if (pmb->packages.AllPackages().count("GRMHD")) { + if (params.Get("cancel_U3_" + bname) && full_grmhd_boundary) { + GRMHD::CancelBoundaryU3(rc.get(), domain, coarse); + } + } + // Always call through to the registered boundary function Flag("Apply "+bname+" boundary: "+btype_name); pkg->KBoundaries[bface](rc, coarse); @@ -349,66 +408,69 @@ void KBoundaries::ApplyBoundary(std::shared_ptr> &rc, IndexD return; } - // Delegate EMF boundaries to the B_CT package - // Only until per-variable boundaries available in Parthenon - // Warning: Even though the EMFs are sync'd separately, - // they still can sneak into "real" boundary exchanges, - // so we can't assume their presence means they are alone - auto& emfpack = rc->PackVariables(std::vector{"B_CT.emf"}); - if (emfpack.GetDim(4) > 0) { - if (params.Get("zero_EMF_" + bname)) { - Flag("ZeroEMF_"+bname); - B_CT::ZeroEMF(rc.get(), domain, emfpack, coarse); - EndFlag(); + // Fixes and special cases for face/edge-centered variabes in B_CT + if (pmb->packages.AllPackages().count("B_CT")) { + // Delegate EMF boundaries to the B_CT package + // Only until per-variable boundaries available in Parthenon + // Warning: Even though the EMFs are sync'd separately, + // they still can sneak into "real" boundary exchanges, + // so we can't assume their presence means they are alone + auto& emfpack = rc->PackVariables(std::vector{"B_CT.emf"}); + if (emfpack.GetDim(4) > 0) { + if (params.Get("zero_EMF_" + bname)) { + Flag("ZeroEMF_"+bname); + B_CT::ZeroBoundaryEMF(rc.get(), domain, emfpack, coarse); + EndFlag(); + } + if (params.Get("average_EMF_" + bname)) { + Flag("AverageEMF_"+bname); + B_CT::AverageBoundaryEMF(rc.get(), domain, emfpack, coarse); + EndFlag(); + } } - if (params.Get("average_EMF_" + bname)) { - Flag("AverageEMF_"+bname); - B_CT::AverageEMF(rc.get(), domain, emfpack, coarse); + + // Correct Parthenon's reflecting conditions on the corresponding face + // Note these are REFLECTING SPECIFIC, not suitable for the similar op w/transmitting + // TODO move this to Parthenon. Move out of this case if we gain non-B_CT face variables + auto fpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost, Metadata::GetUserFlag("SplitVector")}); + if (params.Get("reflect_face_vector_" + bname) && fpack.GetDim(4) > 0) { + Flag("ReflectFace_"+bname); + const TopologicalElement face = FaceOf(bdir); + auto b = KDomain::GetBoundaryRange(rc, domain, face, coarse); + // Zero the last physical face, otherwise invert. + auto i_f = (binner) ? b.ie : b.is; + auto j_f = (binner) ? b.je : b.js; + auto k_f = (binner) ? b.ke : b.ks; + pmb->par_for( + "reflect_face_vector_" + bname, 0, fpack.GetDim(4)-1, b.ks, b.ke, b.js, b.je, b.is, b.ie, + KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) { + const int kk = (bdir == 3) ? k_f - (k - k_f) : k; + const int jj = (bdir == 2) ? j_f - (j - j_f) : j; + const int ii = (bdir == 1) ? i_f - (i - i_f) : i; + fpack(face, v, k, j, i) = ((bdir == 1 && i == i_f) || + (bdir == 2 && j == j_f) || + (bdir == 3 && k == k_f)) ? 0. : -fpack(face, v, kk, jj, ii); + } + ); EndFlag(); } - } - - // Correct Parthenon's reflecting conditions on the corresponding face - // TODO honor SplitVector here, then move it all to Parthenon - auto fpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost}); - if (params.Get("reflect_face_vector_" + bname) && fpack.GetDim(4) > 0) { - Flag("ReflectFace_"+bname); - const TopologicalElement face = FaceOf(bdir); - auto b = KDomain::GetBoundaryRange(rc, domain, face, coarse); - // Zero the last physical face, otherwise invert. - auto i_f = (binner) ? b.ie : b.is; - auto j_f = (binner) ? b.je : b.js; - auto k_f = (binner) ? b.ke : b.ks; - pmb->par_for( - "reflect_face_vector_" + bname, b.ks, b.ke, b.js, b.je, b.is, b.ie, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - const int kk = (bdir == 3) ? k_f - (k - k_f) : k; - const int jj = (bdir == 2) ? j_f - (j - j_f) : j; - const int ii = (bdir == 1) ? i_f - (i - i_f) : i; - fpack(face, 0, k, j, i) = ((bdir == 1 && i == i_f) || - (bdir == 2 && j == j_f) || - (bdir == 3 && k == k_f)) ? 0. : -fpack(face, 0, kk, jj, ii); - } - ); - EndFlag(); - } - // Correct orthogonal B field component to eliminate divergence in last rank - // and ghosts. Used for outflow conditions when field lines will exit domain - if (params.Get("clean_face_B_" + bname) && fpack.GetDim(4) > 0) { - Flag("CleanFaceB_"+bname); - B_CT::DestructiveBoundaryClean(rc.get(), domain, fpack, coarse); - EndFlag(); + // Correct orthogonal B field component to eliminate divergence in last rank + // and ghosts. Used for outflow conditions when field lines will exit domain + auto bfpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost, Metadata::GetUserFlag("B_CT")}); + if (params.Get("clean_face_B_" + bname) && bfpack.GetDim(4) > 0) { + Flag("CleanFaceB_"+bname); + B_CT::DestructiveBoundaryClean(rc.get(), domain, bfpack, coarse); + EndFlag(); + } } - // This function is called in 2 places we might not expect, - // where we still may want to control the physical bounds: + // This function, ApplyBoundary, is called in 2 places we might not expect: // 1. Syncing only the EMF during runs with CT // 2. Syncing boundaries while solving for B field - // but, anything beyond here is really only for the expected case, - // with all the fluid variables - PackIndexMap prims_map; - if (GRMHD::PackMHDPrims(rc.get(), prims_map).GetDim(4) == 0) { + // The above operations are general to these cases + // But, anything beyond this point isn't needed for those cases & may crash + if (!full_grmhd_boundary) { EndFlag(); return; } @@ -535,7 +597,7 @@ TaskStatus KBoundaries::FixFlux(MeshData *md) // valid cell in the face direction. That is, e.g. F1 is valid on // an (N1+1)xN2xN3 grid, F2 on N1x(N2+1)xN3, etc. // These functions do *not* need an extra row outside the domain, - // like B_FluxCT::FixBoundaryFlux does. + // like B_FluxCT::ZeroBoundaryFlux does. const int ndim = pmesh->ndim; // Entire range const IndexRange ibe = pmb0->cellbounds.GetBoundsI(IndexDomain::entire); diff --git a/kharma/boundaries/one_block_transmit.cpp b/kharma/boundaries/one_block_transmit.cpp index dde9492f..1381e540 100644 --- a/kharma/boundaries/one_block_transmit.cpp +++ b/kharma/boundaries/one_block_transmit.cpp @@ -45,20 +45,20 @@ void KBoundaries::TransmitImpl(MeshBlockData *rc, BoundaryFace bface, bool { // Get all cell-centered ghosts, minus anything just used at startup using FC = Metadata::FlagCollection; - FC ghost_vars = FC({Metadata::FillGhost, Metadata::Conserved}) - + FC({Metadata::FillGhost, Metadata::GetUserFlag("Primitive")}) + FC ghost_vars = FC({Metadata::FillGhost, Metadata::Cell}) - FC({Metadata::GetUserFlag("StartupOnly")}); - auto q = rc->PackVariables(ghost_vars, coarse); - TransmitSetTE(rc, q, bface, coarse, false); + PackIndexMap bounds_map; + auto q = rc->PackVariables(ghost_vars, bounds_map, coarse); + TransmitSetTE(rc, q, bface, bounds_map, coarse, false); FC ghost_vars_f = FC({Metadata::FillGhost, Metadata::Face}) - - FC({Metadata::GetUserFlag("StartupOnly")}); - auto q_f = rc->PackVariables(ghost_vars_f, coarse); - TransmitSetTE(rc, q_f, bface, coarse, true); + - FC({Metadata::GetUserFlag("StartupOnly")}); + auto q_f = rc->PackVariables(ghost_vars_f, bounds_map, coarse); + TransmitSetTE(rc, q_f, bface, bounds_map, coarse, true); } -void KBoundaries::TransmitSetTE(MeshBlockData *rc, VariablePack &q, - BoundaryFace bface, bool coarse, bool do_face) +void KBoundaries::TransmitSetTE(MeshBlockData *rc, VariablePack &q, BoundaryFace bface, + PackIndexMap &bounds_map, bool coarse, bool do_face) { // We're sometimes called without any variables to sync (e.g. syncing flags, EMFs), just return if (q.GetDim(4) == 0) return; @@ -100,8 +100,23 @@ void KBoundaries::TransmitSetTE(MeshBlockData *rc, VariablePack &q, // Calculate j just like for reflecting conditions const int jpivot = (binner) ? b.je : b.js; // B3 component on X3 face should be inverted even if not marked "vector" - // TODO honor SplitVector here rather than always inverting + // TODO honor SplitVector and vector components here rather than hard-coding const bool do_face_invert = (el == F3); + // TODO figure out fixing '.vector_component' in Parthenon + const int x3_index_1 = bounds_map["cons.uvec"].second; + const int x3_index_2 = bounds_map["prims.uvec"].second; + const int x3_index_3 = bounds_map["cons.B"].second; + const int x3_index_4 = bounds_map["prims.B"].second; + + //std::cerr << "Transmitting boundary. el=" << static_cast(el) << " do_invert=" << do_face_invert << " corr_face=" << corresponding_face << std::endl; + //std::cerr << "Elements: " << q.GetDim(4) << " vector components: "; + //pmb->par_for( + // "print_" + bname, 0, q.GetDim(4)-1, + // KOKKOS_LAMBDA (const int &v) { + // printf("%d ", q(el, v).vector_component); + // } + //); + //std::cerr << std::endl; pmb->par_for( "transmitting_polar_boundary_" + bname, 0, q.GetDim(4)-1, b.ks, b.ke, b.js, b.je, b.is, b.ie, @@ -109,7 +124,11 @@ void KBoundaries::TransmitSetTE(MeshBlockData *rc, VariablePack &q, const int ki = ((k - ksp + Nk3p2) % Nk3p) + ksp; const int ji = jpivot + reflect_offset + (jpivot - j); const int ii = i; - const Real invert = (do_face_invert || q(el, v).vector_component == X3DIR) ? -1. : 1.; + const Real invert = (do_face_invert || v == x3_index_1 || v == x3_index_2 || + v == x3_index_3 || v == x3_index_4 || + q(el, v).vector_component == X3DIR) ? -1. : 1.; + //if (i == 10 && j == 3 && k == 10) + // printf("Set el %d v %d zone %d %d %d from %d %d %d invert: %f\n", el, v, k, j, i, ki, ji, ii, invert); q(el, v, k, j, i) = (corresponding_face && j == jpivot) ? 0. : invert * q(el, v, ki, ji, ii); } ); diff --git a/kharma/boundaries/one_block_transmit.hpp b/kharma/boundaries/one_block_transmit.hpp index d082105a..3b0de709 100644 --- a/kharma/boundaries/one_block_transmit.hpp +++ b/kharma/boundaries/one_block_transmit.hpp @@ -39,8 +39,8 @@ namespace KBoundaries { // TODO(BSP) privatize probably void TransmitImpl(MeshBlockData *rc, BoundaryFace bface, bool coarse); -void TransmitSetTE(MeshBlockData *rc, VariablePack &q, - BoundaryFace bface, bool coarse, bool do_face); +void TransmitSetTE(MeshBlockData *rc, VariablePack &q, BoundaryFace bface, + PackIndexMap &bounds_map, bool coarse, bool do_face); template inline void OneBlockTransmit(std::shared_ptr> &rc, bool coarse) diff --git a/kharma/coordinates/gr_coordinates.cpp b/kharma/coordinates/gr_coordinates.cpp index 663c1757..a8c8d987 100644 --- a/kharma/coordinates/gr_coordinates.cpp +++ b/kharma/coordinates/gr_coordinates.cpp @@ -96,7 +96,7 @@ GRCoordinates::GRCoordinates(const GRCoordinates &src, int coarsen): UniformCart void init_GRCoordinates(GRCoordinates& G) { const int n1 = G.n1; const int n2 = G.n2; - const int n3 = G.n3; + //const int n3 = G.n3; const bool correct_connections = G.correct_connections; const int connection_average_points = G.connection_average_points; diff --git a/kharma/domain.hpp b/kharma/domain.hpp index a45e775c..2743b000 100644 --- a/kharma/domain.hpp +++ b/kharma/domain.hpp @@ -118,9 +118,9 @@ inline IndexRange3 GetRange(T data, IndexDomain domain, TopologicalElement el=CC const IndexRange ibe = cellbounds.GetBoundsI(IndexDomain::entire, el); const IndexRange jbe = cellbounds.GetBoundsJ(IndexDomain::entire, el); const IndexRange kbe = cellbounds.GetBoundsK(IndexDomain::entire, el); - 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)}; + return IndexRange3{m::max(il.s, ibe.s), m::min(il.e, ibe.e), + m::max(jl.s, jbe.s), m::min(jl.e, jbe.e), + m::max(kl.s, kbe.s), m::min(kl.e, kbe.e)}; } template inline IndexRange3 GetRange(T data, IndexDomain domain, int left_halo, int right_halo, bool coarse=false) @@ -167,23 +167,23 @@ inline IndexRange3 GetPhysicalRange(MeshBlockData* rc) const auto pmb = rc->GetBlockPointer(); return IndexRange3{IsPhysicalBoundary(pmb, BoundaryFace::inner_x1) - ? (uint) bounds.is(IndexDomain::interior) - : (uint) bounds.is(IndexDomain::entire), + ? bounds.is(IndexDomain::interior) + : bounds.is(IndexDomain::entire), IsPhysicalBoundary(pmb, BoundaryFace::outer_x1) - ? (uint) bounds.ie(IndexDomain::interior) - : (uint) bounds.ie(IndexDomain::entire), + ? bounds.ie(IndexDomain::interior) + : bounds.ie(IndexDomain::entire), IsPhysicalBoundary(pmb, BoundaryFace::inner_x2) - ? (uint) bounds.js(IndexDomain::interior) - : (uint) bounds.js(IndexDomain::entire), + ? bounds.js(IndexDomain::interior) + : bounds.js(IndexDomain::entire), IsPhysicalBoundary(pmb, BoundaryFace::outer_x2) - ? (uint) bounds.je(IndexDomain::interior) - : (uint) bounds.je(IndexDomain::entire), + ? bounds.je(IndexDomain::interior) + : bounds.je(IndexDomain::entire), IsPhysicalBoundary(pmb, BoundaryFace::inner_x3) - ? (uint) bounds.ks(IndexDomain::interior) - : (uint) bounds.ks(IndexDomain::entire), + ? bounds.ks(IndexDomain::interior) + : bounds.ks(IndexDomain::entire), IsPhysicalBoundary(pmb, BoundaryFace::outer_x3) - ? (uint) bounds.ke(IndexDomain::interior) - : (uint) bounds.ke(IndexDomain::entire)}; + ? bounds.ke(IndexDomain::interior) + : bounds.ke(IndexDomain::entire)}; } template diff --git a/kharma/driver/imex_step.cpp b/kharma/driver/imex_step.cpp index 8afa1d43..524661b7 100644 --- a/kharma/driver/imex_step.cpp +++ b/kharma/driver/imex_step.cpp @@ -73,6 +73,8 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta const bool use_implicit = pkgs.count("Implicit"); const bool use_jcon = pkgs.count("Current"); const bool use_linesearch = (use_implicit) ? pkgs.at("Implicit")->Param("linesearch") : false; + const bool emhd_enabled = pkgs.count("EMHD"); + const bool use_ideal_guess = (emhd_enabled) ? pkgs.at("GRMHD")->Param("ideal_guess") : false; // Allocate/copy the things we need // TODO these can now be reduced by including the var lists/flags which actually need to be allocated @@ -191,9 +193,20 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta std::vector{Metadata::GetUserFlag("Explicit"), Metadata::Independent}, use_b_ct, stage); + // Update ideal HD variables so that they can be used as a guess for the solver. + // Only used if emhd/ideal_guess is enabled. + // An additional `AddStateUpdate` task just for variables marked with the `IdealGuess` flag + auto t_ideal_guess = t_update; + if (use_ideal_guess) { + t_ideal_guess = KHARMADriver::AddStateUpdateIdealGuess(t_sources, tl, md_full_step_init.get(), md_sub_step_init.get(), + md_flux_src.get(), md_solver.get(), + std::vector{Metadata::GetUserFlag("IdealGuess"), Metadata::Independent}, + false, stage); + } + // Make sure the primitive values of *explicitly-evolved* variables are updated. // Packages with implicitly-evolved vars should only register BoundaryUtoP or BoundaryPtoU - auto t_explicit_UtoP = tl.AddTask(t_update, Packages::MeshUtoP, md_solver.get(), IndexDomain::entire, false); + auto t_explicit_UtoP = tl.AddTask(t_ideal_guess, Packages::MeshUtoP, md_solver.get(), IndexDomain::entire, false); // Done with explicit update auto t_explicit = t_explicit_UtoP; @@ -204,9 +217,17 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta std::shared_ptr> &md_linesearch = (use_linesearch) ? pmesh->mesh_data.GetOrAdd("linesearch", i) : md_solver; // Copy the current state of any implicitly-evolved vars (at least the prims) in as a guess. - // This sets md_solver = md_sub_step_init - auto t_copy_guess = tl.AddTask(t_sources, Copy>, std::vector({Metadata::GetUserFlag("Implicit")}), + // If we aren't using the ideal solution as the guess, set md_solver = md_sub_step_init + auto t_copy_guess = t_explicit; + auto t_copy_emhd_vars = t_explicit; + if (use_ideal_guess) { + t_copy_emhd_vars = tl.AddTask(t_sources, Copy>, std::vector({Metadata::GetUserFlag("EMHDVar")}), md_sub_step_init.get(), md_solver.get()); + t_copy_guess = t_copy_emhd_vars; + } else { + t_copy_guess = tl.AddTask(t_sources, Copy>, std::vector({Metadata::GetUserFlag("Implicit")}), + md_sub_step_init.get(), md_solver.get()); + } auto t_guess_ready = t_explicit | t_copy_guess; @@ -219,7 +240,6 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta md_solver.get(), md_linesearch.get()); } - // Time-step implicit variables by root-finding the residual. // This calculates the primitive values after the substep for all "isImplicit" variables -- // no need for separately adding the flux divergence or calling UtoP diff --git a/kharma/driver/kharma_driver.cpp b/kharma/driver/kharma_driver.cpp index c187d72b..1878d760 100644 --- a/kharma/driver/kharma_driver.cpp +++ b/kharma/driver/kharma_driver.cpp @@ -82,6 +82,14 @@ std::shared_ptr KHARMADriver::Initialize(ParameterInput *pin, std // so we define the flags here to avoid loading order issues Metadata::AddUserFlag("Implicit"); Metadata::AddUserFlag("Explicit"); + // Add a flag if we wish to use ideal variables explicitly evolved as guess for implicit update. + // GRIM uses the fluid state of the previous (sub-)step. + // The logic here is that the non-ideal variables do not significantly contribute to the + // stress-energy tensor. We can explicitly update the ideal MHD variables and hope that this + // takes us to a region in the parameter space of state variables that + // is close to the true solution. The corrections obtained from the implicit update would then + // be small corrections. + Metadata::AddUserFlag("IdealGuess"); // 1. One flag to mark the primitive variables specifically // (Parthenon has Metadata::Conserved already) @@ -290,7 +298,8 @@ TaskID KHARMADriver::AddFOFC(TaskID& t_start, TaskList& tl, MeshData *md, auto pmb0 = md->GetBlockData(0)->GetBlockPointer(); auto& pkgs = pmb0->packages.AllPackages(); - const Floors::Prescription fofc_floors = pmb0->packages.Get("Flux")->Param("fofc_prescription"); + const Floors::Prescription fofc_floors = pmb0->packages.Get("Flux")->Param("fofc_prescription"); + const Floors::Prescription fofc_floors_inner = pmb0->packages.Get("Flux")->Param("fofc_prescription_inner"); // Populate guess source term with divergence of the existing fluxes // NOTE this does not include source terms! Though, could call them here tbh @@ -313,7 +322,7 @@ TaskID KHARMADriver::AddFOFC(TaskID& t_start, TaskList& tl, MeshData *md, auto t_guess_Bp = tl.AddTask(t_guess_update, B_FluxCT::MeshUtoP, guess, IndexDomain::entire, false); auto t_guess_prims = tl.AddTask(t_guess_Bp, Inverter::MeshUtoP, guess, IndexDomain::entire, false); // Check and mark floors - auto t_mark_floors = tl.AddTask(t_guess_prims, Floors::DetermineGRMHDFloors, guess, IndexDomain::entire, fofc_floors); + auto t_mark_floors = tl.AddTask(t_guess_prims, Floors::DetermineGRMHDFloors, guess, IndexDomain::entire, fofc_floors, fofc_floors_inner); // Determine which cells are FOFC in our block auto t_mark_fofc = tl.AddTask(t_mark_floors, Flux::MarkFOFC, guess); // Sync with neighbor blocks. This seems to ameliorate an increasing divB on X1 boundaries in GR, @@ -377,6 +386,56 @@ TaskID KHARMADriver::AddStateUpdate(TaskID& t_start, TaskList& tl, MeshData *md_full_step_init, + MeshData *md_sub_step_init, MeshData *md_flux_src, MeshData *md_update, + std::vector flags, bool update_face, int stage) +{ + // Update variables marked as IdealGuess + // Add any proportion of the step start required by the integrator (e.g., RK2) + std::vector flags_cell = flags; flags_cell.push_back(Metadata::Cell); + std::vector flags_face = flags; flags_face.push_back(Metadata::Face); + // TODO splitting this is stupid, but maybe the parallelization actually helps? Eh. + auto t_avg_data_c = tl.AddTask(t_start, Update::WeightedSumData, MeshData>, + std::vector(flags_cell), + md_sub_step_init, md_full_step_init, + integrator->gam0[stage-1], integrator->gam1[stage-1], + md_update); + auto t_avg_data = t_avg_data_c; + if (update_face) { + t_avg_data = tl.AddTask(t_start, WeightedSumDataFace, + std::vector(flags_face), + md_sub_step_init, md_full_step_init, + integrator->gam0[stage-1], integrator->gam1[stage-1], + md_update); + } + // apply du/dt to the result + auto t_update_c = tl.AddTask(t_avg_data, Update::WeightedSumData, MeshData>, + std::vector(flags_cell), + md_update, md_flux_src, + 1.0, integrator->beta[stage-1] * integrator->dt, + md_update); + auto t_update = t_update_c; + if (update_face) { + t_update = tl.AddTask(t_avg_data, WeightedSumDataFace, + std::vector(flags_face), + md_update, md_flux_src, + 1.0, integrator->beta[stage-1] * integrator->dt, + md_update); + } + + // We'll be running UtoP after this, which needs a guess in order to converge, so we copy in md_sub_step_init + auto t_copy_prims = t_update; + auto pmb0 = md_full_step_init->GetBlockData(0)->GetBlockPointer(); + auto& pkgs = pmb0->packages.AllPackages(); + if (pkgs.at("GRMHD")->Param("ideal_guess")) { + t_copy_prims = tl.AddTask(t_start, Copy>, + std::vector({Metadata::GetUserFlag("HD"), Metadata::GetUserFlag("Primitive")}), + md_sub_step_init, md_update); + } + + return t_copy_prims | t_update; +} + void KHARMADriver::SetGlobalTimeStep() { // TODO(BSP) apply the limits from GRMHD package here diff --git a/kharma/driver/kharma_driver.hpp b/kharma/driver/kharma_driver.hpp index d660e245..42abf273 100644 --- a/kharma/driver/kharma_driver.hpp +++ b/kharma/driver/kharma_driver.hpp @@ -133,6 +133,15 @@ class KHARMADriver : public MultiStageDriver { MeshData *md_flux_src, MeshData *md_update, std::vector flags, bool update_face, int stage); + /** + * This function updates a state md_update with the results of an explicit source term calculation + * placed in md_flux_src. It is similar to `AddStateUpdate` but applies only to variables marked + * with the `IdealGuess` flag. + */ + TaskID AddStateUpdateIdealGuess(TaskID& t_start, TaskList& tl, MeshData *md_full_step_init, MeshData *md_sub_step_init, + MeshData *md_flux_src, MeshData *md_update, std::vector flags, + bool update_face, int stage); + /** * Add a synchronization retion to an existing TaskCollection tc. * Since the region is self-contained, does not return a TaskID diff --git a/kharma/electrons/electrons.cpp b/kharma/electrons/electrons.cpp index 3b1601a6..eeecd134 100644 --- a/kharma/electrons/electrons.cpp +++ b/kharma/electrons/electrons.cpp @@ -254,8 +254,6 @@ void BlockUtoP(MeshBlockData *rc, IndexDomain domain, bool coarse) // And then the local density GridScalar rho_U = rc->Get("cons.rho").data; - const auto& G = pmb->coords; - auto bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; int is = bounds.is(domain), ie = bounds.ie(domain); int js = bounds.js(domain), je = bounds.je(domain); @@ -563,7 +561,8 @@ void ApplyFloors(MeshBlockData *mbd, IndexDomain domain) const auto& G = pmb->coords; const Real gam = packages.Get("GRMHD")->Param("gamma"); - const Floors::Prescription floors = packages.Get("Floors")->Param("prescription"); + const Floors::Prescription floors = packages.Get("Floors")->Param("prescription"); + const Floors::Prescription floors_inner = packages.Get("Floors")->Param("prescription_inner"); const IndexRange3 b = KDomain::GetRange(mbd, domain); pmb->par_for("apply_electrons_floors", b.ks, b.ke, b.js, b.je, b.is, b.ie, @@ -571,9 +570,19 @@ void ApplyFloors(MeshBlockData *mbd, IndexDomain domain) // Also apply the ceiling to the advected entropy KTOT, if we're keeping track of that // (either for electrons, or robust primitive inversions in future) - if (m_p.KTOT >= 0 && (P(m_p.KTOT, k, j, i) > floors.ktot_max)) { - fflag(0, k, j, i) = Floors::FFlag::KTOT | (int) fflag(0, k, j, i); - P(m_p.KTOT, k, j, i) = floors.ktot_max; + Real ktot_max; + if (m_p.KTOT >= 0) { + if (floors.radius_dependent_floors && G.coords.is_spherical() + && G.r(k, j, i) < floors.floors_switch_r) { + ktot_max = floors_inner.ktot_max; + } 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; + } } // TODO(BSP) restore Ressler adjustment option diff --git a/kharma/floors/floors.cpp b/kharma/floors/floors.cpp index e0aa07d5..c0685421 100644 --- a/kharma/floors/floors.cpp +++ b/kharma/floors/floors.cpp @@ -89,10 +89,16 @@ std::shared_ptr Floors::Initialize(ParameterInput *pin, std::shar params.Add("frame_switch_beta", frame_switch_beta); } - // Disable all floors. It is obviously tremendously inadvisable to do this - // Note this will only be recorded if false, otherwise we're not loaded at all! - bool disable_floors = pin->GetOrAddBoolean("floors", "disable_floors", false); - params.Add("disable_floors", disable_floors); + // Radius-dependent floors and ceiling in a given frame + // Default presciption struct refers to outer domain, i.e., beyond 'floor_switch_r' + // Make an additional prescription struct for inner domain. + // There are two Floors::Prescription objects even if there is no radius-dependence, + // the values will simply be the same if radius-dependent floors are not enabled. + // Avoids a bunch of if (radius_dependent_floors) else while determining floors. + if (pin->DoesBlockExist("floors_inner")) + params.Add("prescription_inner", MakePrescriptionInner(pin, MakePrescription(pin))); + else + params.Add("prescription_inner", MakePrescriptionInner(pin, MakePrescription(pin)), "floors"); // These preserve floor values between the "mark" pass and the actual floor application // We need them even if floors are disabled, to apply initial values based on some prescription @@ -149,6 +155,8 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData * // JUST rho & u geometric floors_tmp.rho_min_geom = pin->GetOrAddReal("floors", "rho_min_geom", 1e-6); floors_tmp.u_min_geom = pin->GetOrAddReal("floors", "u_min_geom", 1e-8); + floors_tmp.rho_min_const = pin->GetOrAddReal("floors", "rho_min_const", 1e-20); + floors_tmp.u_min_const = pin->GetOrAddReal("floors", "u_min_const", 1e-20); // Disable everything else, even if it's specified floors_tmp.bsq_over_rho_max = 1e20; @@ -170,10 +178,11 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData * pmb->par_for("apply_initial_floors", b.ks, b.ke, b.js, b.je, b.is, b.ie, KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { Real rhoflr_max, uflr_max; - int fflag = determine_floors(G, P, m_p, gam, k, j, i, floors, rhoflr_max, uflr_max); + // Initial floors, so the radius-dependence of floors don't matter that much. + int fflag = determine_floors(G, P, m_p, gam, k, j, i, floors, floors, rhoflr_max, uflr_max); if (fflag) { apply_floors(G, P, m_p, gam, k, j, i, rhoflr_max, uflr_max, U, m_u); - apply_ceilings(G, P, m_p, gam, k, j, i, floors, U, m_u); + apply_ceilings(G, P, m_p, gam, k, j, i, floors, floors, U, m_u); // P->U for any modified zones Flux::p_to_u_mhd(G, P, m_p, emhd_params, gam, k, j, i, U, m_u, Loci::center); } @@ -184,7 +193,8 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData * return TaskStatus::complete; } -TaskStatus Floors::DetermineGRMHDFloors(MeshData *md, IndexDomain domain, const Floors::Prescription& floors) +TaskStatus Floors::DetermineGRMHDFloors(MeshData *md, IndexDomain domain, + const Floors::Prescription& floors, const Floors::Prescription& floors_inner) { auto pmb0 = md->GetBlockData(0)->GetBlockPointer(); @@ -208,7 +218,7 @@ TaskStatus Floors::DetermineGRMHDFloors(MeshData *md, IndexDomain domain, KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i) { const auto& G = P.GetCoords(b); fflag(b, 0, k, j, i) = static_cast(fflag(b, 0, k, j, i)) | - determine_floors(G, P(b), m_p, gam, k, j, i, floors, + determine_floors(G, P(b), m_p, gam, k, j, i, floors, floors_inner, floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i)); } ); diff --git a/kharma/floors/floors.hpp b/kharma/floors/floors.hpp index 9531281c..6e2b2363 100644 --- a/kharma/floors/floors.hpp +++ b/kharma/floors/floors.hpp @@ -111,6 +111,9 @@ class Prescription { Real gamma_max; // Floor options (frame was MOVED to templating) bool use_r_char, temp_adjust_u, adjust_k; + // Radius dependent floors? + bool radius_dependent_floors; + Real floors_switch_r; }; inline Prescription MakePrescription(parthenon::ParameterInput *pin, std::string block="floors") @@ -157,9 +160,60 @@ inline Prescription MakePrescription(parthenon::ParameterInput *pin, std::string // Limit the fluid Lorentz factor gamma p.gamma_max = pin->GetOrAddReal(block, "gamma_max", 50.); + p.radius_dependent_floors = pin->GetOrAddBoolean("floors", "radius_dependent_floors", false); + p.floors_switch_r = pin->GetOrAddReal("floors", "floors_switch_r", 50.); + return p; } +/** + * Set prescription struct for inner domain (r < floor_switch_r) if 'radius_dependent_floors' is enabled. + * Sets values provided in the 'floors_inner' block in input par file if provided, + * else sets it equal to the values in the whole/outer domain. + */ +inline Prescription MakePrescriptionInner(parthenon::ParameterInput *pin, Prescription p_outer, std::string block="floors_inner") +{ + // TODO(BSP) I wonder if there's an easier way to "set if parameter exists" from pin, that would be broadly useful + Prescription p_inner; + + // Floor parameters + if (pin->GetBoolean("coordinates", "spherical")) { + // In spherical systems, floors drop as r^2, so set them higher by default + p_inner.rho_min_geom = pin->GetOrAddReal(block, "rho_min_geom", p_outer.rho_min_geom); + p_inner.u_min_geom = pin->GetOrAddReal(block, "u_min_geom", p_outer.u_min_geom); + // Some constant for large distances. New, out of the way by default + p_inner.rho_min_const = pin->GetOrAddReal(block, "rho_min_const", p_outer.rho_min_const); + p_inner.u_min_const = pin->GetOrAddReal(block, "u_min_const", p_outer.u_min_const); + } else { // TODO spherical cart will also have both + // Accept old names + Real rho_min_const_default = pin->DoesParameterExist(block, "rho_min_geom") ? + pin->GetReal(block, "rho_min_geom") : p_outer.rho_min_geom; + Real u_min_const_default = pin->DoesParameterExist(block, "u_min_geom") ? + pin->GetReal(block, "u_min_geom") : p_outer.u_min_geom; + p_inner.rho_min_const = pin->GetOrAddReal(block, "rho_min_const", rho_min_const_default); + p_inner.u_min_const = pin->GetOrAddReal(block, "u_min_const", u_min_const_default); + } + + p_inner.use_r_char = pin->GetOrAddBoolean(block, "use_r_char", p_outer.use_r_char); + p_inner.r_char = pin->GetOrAddReal(block, "r_char", p_outer.r_char); + + p_inner.bsq_over_rho_max = pin->GetOrAddReal(block, "bsq_over_rho_max", p_outer.bsq_over_rho_max); + p_inner.bsq_over_u_max = pin->GetOrAddReal(block, "bsq_over_u_max", p_outer.bsq_over_u_max); + + p_inner.u_over_rho_max = pin->GetOrAddReal(block, "u_over_rho_max", p_outer.u_over_rho_max); + p_inner.ktot_max = pin->GetOrAddReal(block, "ktot_max", p_outer.ktot_max); + p_inner.temp_adjust_u = pin->GetOrAddBoolean(block, "temp_adjust_u", p_outer.temp_adjust_u); + p_inner.adjust_k = pin->GetOrAddBoolean(block, "adjust_k", p_outer.adjust_k); + + p_inner.gamma_max = pin->GetOrAddReal(block, "gamma_max", p_outer.gamma_max); + + // Always grab these from p_outer, they should never differ between outer/inner floors + p_inner.radius_dependent_floors = p_outer.radius_dependent_floors; + p_inner.floors_switch_r = p_outer.floors_switch_r; + + return p_inner; +} + /** * Initialization. Set parameters. */ @@ -168,8 +222,8 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

U */ @@ -181,7 +235,8 @@ TaskStatus ApplyGRMHDFloors(MeshData *md, IndexDomain domain); * 2. fflag, which floors were hit by the current state * This is what ApplyFloors uses to determine the floor values/locations */ -TaskStatus DetermineGRMHDFloors(MeshData *md, IndexDomain domain, const Floors::Prescription& floors); +TaskStatus DetermineGRMHDFloors(MeshData *md, IndexDomain domain, + const Floors::Prescription& floors, const Floors::Prescription& floors_inner); /** * Apply the same floors as above, in the same way, except: @@ -189,7 +244,7 @@ TaskStatus DetermineGRMHDFloors(MeshData *md, IndexDomain domain, const Fl * 2. Don't record results to 'fflag' or 'pflag' * Used for problems where some part of the domain is initialized to * "whatever the floor value is." - * *This function can be called even if the Floors package is not initialized.* + * *This function can be called even if the Floors package is not initialized, and ignores "floors/on=false"* */ TaskStatus ApplyInitialFloors(ParameterInput *pin, MeshBlockData *mbd, IndexDomain domain); diff --git a/kharma/floors/floors_functions.hpp b/kharma/floors/floors_functions.hpp index 090e6087..8a9fe905 100644 --- a/kharma/floors/floors_functions.hpp +++ b/kharma/floors/floors_functions.hpp @@ -43,18 +43,27 @@ namespace Floors { /** * Apply all ceilings together, currently at most one on velocity and two on internal energy + * TODO REALLY need to take fflag here and only compute existing values if we know the floor was hit * * LOCKSTEP: this function respects P and returns consistent P<->U */ KOKKOS_INLINE_FUNCTION void apply_ceilings(const GRCoordinates& G, const VariablePack& P, const VarMap& m_p, - const Real& gam, const int& k, const int& j, const int& i, const Floors::Prescription& floors, + const Real& gam, const int& k, const int& j, const int& i, + const Floors::Prescription& floors, const Floors::Prescription& floors_inner, const VariablePack& U, const VarMap& m_u, const Loci loc=Loci::center) { - // First apply ceilings: + // Choose our floor scheme + const Floors::Prescription& myfloors = (floors.radius_dependent_floors && G.coords.is_spherical() + && G.r(k, j, i) < floors.floors_switch_r) ? floors_inner : floors; + + // Compute max values for ceilings + Real gamma = GRMHD::lorentz_calc(G, P, m_p, k, j, i, loc); + Real ktot = (gam - 1.) * P(m_p.UU, k, j, i) / m::pow(P(m_p.RHO, k, j, i), gam); + Real u_over_rho = P(m_p.UU, k, j, i) / P(m_p.RHO, k, j, i); + // 1. Limit gamma with respect to normal observer - Real gamma = GRMHD::lorentz_calc(G, P, m_p, k, j, i, loc); - if (gamma > floors.gamma_max) { - Real f = m::sqrt((SQR(floors.gamma_max) - 1.) / (SQR(gamma) - 1.)); + if (gamma > myfloors.gamma_max) { + Real f = m::sqrt((SQR(myfloors.gamma_max) - 1.) / (SQR(gamma) - 1.)); VLOOP P(m_p.U1+v, k, j, i) *= f; } @@ -62,21 +71,24 @@ KOKKOS_INLINE_FUNCTION void apply_ceilings(const GRCoordinates& G, const Variabl // Note this technically applies the condition *one step sooner* than legacy, since it operates on // the entropy as calculated from current conditions, rather than the value kept from the previous // step for calculating dissipation. - Real ktot = (gam - 1.) * P(m_p.UU, k, j, i) / m::pow(P(m_p.RHO, k, j, i), gam); - if (ktot > floors.ktot_max) { - P(m_p.UU, k, j, i) = floors.ktot_max / ktot * P(m_p.UU, k, j, i); + if (ktot > myfloors.ktot_max) { + P(m_p.UU, k, j, i) = myfloors.ktot_max / ktot * P(m_p.UU, k, j, i); } // 3. Limit the temperature by controlling u. Can optionally add density instead, implemented in apply_floors - if (floors.temp_adjust_u && P(m_p.UU, k, j, i) / P(m_p.RHO, k, j, i) > floors.u_over_rho_max) { - P(m_p.UU, k, j, i) = floors.u_over_rho_max * P(m_p.RHO, k, j, i); + if (myfloors.temp_adjust_u && P(m_p.UU, k, j, i) / P(m_p.RHO, k, j, i) > myfloors.u_over_rho_max) { + P(m_p.UU, k, j, i) = myfloors.u_over_rho_max * P(m_p.RHO, k, j, i); } } KOKKOS_INLINE_FUNCTION int determine_floors(const GRCoordinates& G, const VariablePack& P, const VarMap& m_p, const Real& gam, const int& k, const int& j, const int& i, const Floors::Prescription& floors, - Real& rhoflr_max, Real& uflr_max) + const Floors::Prescription& floors_inner, Real& rhoflr_max, Real& uflr_max) { + // Choose our floor scheme + const Floors::Prescription& myfloors = (floors.radius_dependent_floors && G.coords.is_spherical() + && G.r(k, j, i) < floors.floors_switch_r) ? floors_inner : floors; + // Calculate the different floor values in play: // 1. Geometric hard floors, not based on fluid relationships // TODO(BSP) can this be cached if it's slow? @@ -84,20 +96,23 @@ KOKKOS_INLINE_FUNCTION int determine_floors(const GRCoordinates& G, const Variab if(G.coords.is_spherical()) { const GReal r = G.r(k, j, i); // r_char sets more aggressive floor close to EH but backs off - Real rhoscal = (floors.use_r_char) ? 1. / ((r*r) * (1 + r / floors.r_char)) : 1. / m::sqrt(r*r*r); - rhoflr_geom = m::max(floors.rho_min_geom * rhoscal, floors.rho_min_const); - uflr_geom = m::max(floors.u_min_geom * m::pow(rhoscal, gam), floors.u_min_const); + Real rhoscal = (myfloors.use_r_char) ? 1. / ((r*r) * (1 + r / myfloors.r_char)) : 1. / m::sqrt(r*r*r); + rhoflr_geom = m::max(myfloors.rho_min_geom * rhoscal, myfloors.rho_min_const); + uflr_geom = m::max(myfloors.u_min_geom * m::pow(rhoscal, gam), myfloors.u_min_const); } else { - rhoflr_geom = floors.rho_min_const; - uflr_geom = floors.u_min_const; + rhoflr_geom = myfloors.rho_min_const; + uflr_geom = myfloors.u_min_const; } // 2. Magnetization ceilings: impose maximum magnetization sigma = bsq/rho, and inverse beta prop. to bsq/U FourVectors Dtmp; GRMHD::calc_4vecs(G, P, m_p, k, j, i, Loci::center, Dtmp); - Real bsq = dot(Dtmp.bcon, Dtmp.bcov); - Real rhoflr_b = bsq / floors.bsq_over_rho_max; - Real uflr_b = bsq / floors.bsq_over_u_max; + Real rhoflr_b, uflr_b; + // Radius-dependent floors. + // Used with spherical coordinate system, i.e., G.r(k,j,i) exist + Real bsq = dot(Dtmp.bcon, Dtmp.bcov); + rhoflr_b = bsq / myfloors.bsq_over_rho_max; + uflr_b = bsq / myfloors.bsq_over_u_max; // Evaluate max U floor, needed for temp ceiling below uflr_max = m::max(uflr_geom, uflr_b); @@ -106,10 +121,11 @@ KOKKOS_INLINE_FUNCTION int determine_floors(const GRCoordinates& G, const Variab const auto& u = P(m_p.UU, k, j, i); int fflag = 0; - if (!floors.temp_adjust_u) { + if (!myfloors.temp_adjust_u) { // 3. Temperature ceiling: impose maximum temperature u/rho // Take floors on U into account - double rhoflr_temp = m::max(u, uflr_max) / floors.u_over_rho_max; + const double rhoflr_temp = m::max(u, uflr_max) / myfloors.u_over_rho_max; + // Record hitting temperature ceiling fflag |= (rhoflr_temp > rho) * FFlag::TEMP; @@ -130,14 +146,14 @@ KOKKOS_INLINE_FUNCTION int determine_floors(const GRCoordinates& G, const Variab fflag |= (uflr_b > u) * FFlag::B_U; } - // Then ceilings, need to record these for FOFC. See real implementation. - if (GRMHD::lorentz_calc(G, P, m_p, k, j, i, Loci::center) > floors.gamma_max) + // Then ceilings, need to record these for FOFC. See real implementation for details + if (GRMHD::lorentz_calc(G, P, m_p, k, j, i, Loci::center) > myfloors.gamma_max) fflag |= FFlag::GAMMA; - if ((gam - 1.) * P(m_p.UU, k, j, i) / m::pow(P(m_p.RHO, k, j, i), gam) > floors.ktot_max) + if ((gam - 1.) * P(m_p.UU, k, j, i) / m::pow(P(m_p.RHO, k, j, i), gam) > myfloors.ktot_max) fflag |= FFlag::KTOT; - if (floors.temp_adjust_u && (P(m_p.UU, k, j, i) / P(m_p.RHO, k, j, i) > floors.u_over_rho_max)) + if (myfloors.temp_adjust_u && (P(m_p.UU, k, j, i) / P(m_p.RHO, k, j, i) > myfloors.u_over_rho_max)) fflag |= FFlag::TEMP; return fflag; @@ -310,29 +326,24 @@ KOKKOS_INLINE_FUNCTION int apply_floors(FLOOR template KOKKOS_INLINE_FUNCTION int apply_geo_floors(const GRCoordinates& G, Local& P, const VarMap& m, const Real& gam, const int& j, const int& i, - const Floors::Prescription& floors, const Loci loc=Loci::center) + const Floors::Prescription& floors, const Floors::Prescription& floors_inner, + const Loci loc=Loci::center) { + // Choose our floor scheme + const Floors::Prescription& myfloors = (floors.radius_dependent_floors && G.coords.is_spherical() + && G.r(0, j, i) < floors.floors_switch_r) ? floors_inner : floors; + // Apply only the geometric floors Real rhoflr_geom, uflr_geom; if(G.coords.is_spherical()) { - GReal Xembed[GR_DIM]; - G.coord_embed(0, j, i, loc, Xembed); - GReal r = Xembed[1]; - - if (floors.use_r_char) { - // Steeper floor from iharm3d - Real rhoscal = 1. / ((r*r) * (1 + r / floors.r_char)); - rhoflr_geom = floors.rho_min_geom * rhoscal; - uflr_geom = floors.u_min_geom * m::pow(rhoscal, gam); - } else { - // Original floors from iharm2d - Real rhoscal = 1. / m::sqrt(r*r*r); - rhoflr_geom = floors.rho_min_geom * rhoscal; - uflr_geom = floors.u_min_geom * rhoscal / r; - } + const GReal r = G.r(0, j, i); + // r_char sets more aggressive floor close to EH but backs off + Real rhoscal = (myfloors.use_r_char) ? 1. / ((r*r) * (1 + r / myfloors.r_char)) : 1. / m::sqrt(r*r*r); + rhoflr_geom = m::max(myfloors.rho_min_geom * rhoscal, myfloors.rho_min_const); + uflr_geom = m::max(myfloors.u_min_geom * m::pow(rhoscal, gam), myfloors.u_min_const); } else { - rhoflr_geom = floors.rho_min_geom; - uflr_geom = floors.u_min_geom; + rhoflr_geom = myfloors.rho_min_const; + uflr_geom = myfloors.u_min_const; } int fflag = 0; @@ -349,29 +360,24 @@ KOKKOS_INLINE_FUNCTION int apply_geo_floors(const GRCoordinates& G, Local& P, co template KOKKOS_INLINE_FUNCTION int apply_geo_floors(const GRCoordinates& G, Global& P, const VarMap& m, const Real& gam, const int& k, const int& j, const int& i, - const Floors::Prescription& floors, const Loci loc=Loci::center) + const Floors::Prescription& floors, const Floors::Prescription& floors_inner, + const Loci loc=Loci::center) { + // Choose our floor scheme + const Floors::Prescription& myfloors = (floors.radius_dependent_floors && G.coords.is_spherical() + && G.r(0, j, i) < floors.floors_switch_r) ? floors_inner : floors; + // Apply only the geometric floors Real rhoflr_geom, uflr_geom; if(G.coords.is_spherical()) { - GReal Xembed[GR_DIM]; - G.coord_embed(k, j, i, loc, Xembed); - GReal r = Xembed[1]; - - if (floors.use_r_char) { - // Steeper floor from iharm3d - Real rhoscal = 1. / ((r*r) * (1 + r / floors.r_char)); - rhoflr_geom = floors.rho_min_geom * rhoscal; - uflr_geom = floors.u_min_geom * m::pow(rhoscal, gam); - } else { - // Original floors from iharm2d - Real rhoscal = 1. / m::sqrt(r*r*r); - rhoflr_geom = floors.rho_min_geom * rhoscal; - uflr_geom = floors.u_min_geom * rhoscal / r; - } + const GReal r = G.r(0, j, i); + // r_char sets more aggressive floor close to EH but backs off + Real rhoscal = (myfloors.use_r_char) ? 1. / ((r*r) * (1 + r / myfloors.r_char)) : 1. / m::sqrt(r*r*r); + rhoflr_geom = m::max(myfloors.rho_min_geom * rhoscal, myfloors.rho_min_const); + uflr_geom = m::max(myfloors.u_min_geom * m::pow(rhoscal, gam), myfloors.u_min_const); } else { - rhoflr_geom = floors.rho_min_geom; - uflr_geom = floors.u_min_geom; + rhoflr_geom = myfloors.rho_min_const; + uflr_geom = myfloors.u_min_const; } int fflag = 0; diff --git a/kharma/floors/floors_impl.hpp b/kharma/floors/floors_impl.hpp index 60787b5c..58738127 100644 --- a/kharma/floors/floors_impl.hpp +++ b/kharma/floors/floors_impl.hpp @@ -63,9 +63,10 @@ TaskStatus ApplyFloorsInFrame(MeshData *md, IndexDomain domain) const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb0->packages); // Still needed for ceilings and determining floors const Floors::Prescription floors = pmb0->packages.Get("Floors")->Param("prescription"); + const Floors::Prescription floors_inner = pmb0->packages.Get("Floors")->Param("prescription_inner"); // Determine floors - DetermineGRMHDFloors(md, domain, floors); + DetermineGRMHDFloors(md, domain, floors, floors_inner); const IndexRange3 b = KDomain::GetRange(md, domain); const IndexRange block = IndexRange{0, P.GetDim(5) - 1}; @@ -83,7 +84,7 @@ TaskStatus ApplyFloorsInFrame(MeshData *md, IndexDomain domain) if (pflag_l) pflag(b, 0, k, j, i) = pflag_l; // Apply ceilings *after* floors, to make the temperature ceiling better-behaved - apply_ceilings(G, P(b), m_p, gam, k, j, i, floors, U(b), m_u); + apply_ceilings(G, P(b), m_p, gam, k, j, i, floors, floors_inner, U(b), m_u); // P->U for any modified zones Flux::p_to_u_mhd(G, P(b), m_p, emhd_params, gam, k, j, i, U(b), m_u, Loci::center); diff --git a/kharma/flux/flux.cpp b/kharma/flux/flux.cpp index c9b37d12..76aeba52 100644 --- a/kharma/flux/flux.cpp +++ b/kharma/flux/flux.cpp @@ -197,6 +197,11 @@ std::shared_ptr Flux::Initialize(ParameterInput *pin, std::shared bool use_source_term = pin->GetOrAddBoolean("fofc", "use_source_term", false); params.Add("fofc_use_source_term", use_source_term); + int fofc_polar_cells = pin->GetOrAddInteger("fofc", "polar_cells", 0); + params.Add("fofc_polar_cells", fofc_polar_cells); + const GReal eh_buffer = pin->GetOrAddReal("fofc", "eh_buffer", 0.1); + params.Add("fofc_eh_buffer", eh_buffer); + if (packages->AllPackages().count("B_CT")) { // Use consistent B for FOFC (see above) // It is mildly inadvisable to disable this @@ -208,8 +213,14 @@ std::shared_ptr Flux::Initialize(ParameterInput *pin, std::shared // TODO even post-reconstruction/reconstruction fallback? if (!pin->DoesBlockExist("fofc_floors")) { params.Add("fofc_prescription", Floors::MakePrescription(pin, "floors")); + if (pin->DoesBlockExist("floors_inner")) + params.Add("fofc_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "floors"), "floors_inner")); + else + params.Add("fofc_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "floors"), "floors")); } else { + // Override inner and outer floors with `fofc_floors` block params.Add("fofc_prescription", Floors::MakePrescription(pin, "fofc_floors")); + params.Add("fofc_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "fofc_floors"), "fofc_floors")); } // Flag for whether FOFC was applied, for diagnostics diff --git a/kharma/flux/flux.hpp b/kharma/flux/flux.hpp index 96ea78e3..970b363f 100644 --- a/kharma/flux/flux.hpp +++ b/kharma/flux/flux.hpp @@ -91,6 +91,8 @@ TaskStatus MeshPtoU(MeshData *md, IndexDomain domain, bool coarse=false); /** * As above, except that IndexDomains of ghost cells are taken to cover * cells *sent* (that is, part of the domain) rather than received + * + * UNUSED */ TaskStatus BlockPtoU_Send(MeshBlockData *rc, IndexDomain domain, bool coarse); diff --git a/kharma/flux/fofc.cpp b/kharma/flux/fofc.cpp index 60fda042..e366ddbd 100644 --- a/kharma/flux/fofc.cpp +++ b/kharma/flux/fofc.cpp @@ -54,8 +54,11 @@ TaskStatus Flux::MarkFOFC(MeshData *guess) auto fofcflag = guess->PackVariables(std::vector{"fofcflag"}); // Parameters + const auto& pars = pmb0->packages.Get("Flux")->AllParams(); const bool spherical = pmb0->coords.coords.is_spherical(); const GReal r_eh = pmb0->coords.coords.get_horizon(); + const int polar_cells = pars.Get("fofc_polar_cells"); + const GReal eh_buffer = pars.Get("fofc_eh_buffer"); // Pre-mark cells which will need fluxes reduced. // This avoids a race condition marking them multiple times when iterating faces, @@ -68,13 +71,46 @@ TaskStatus Flux::MarkFOFC(MeshData *guess) // if cell failed to invert or would call floors... // TODO preserve cause in the fofcflag if (static_cast(fflag(b, 0, k, j, i)) || //Inverter::failed(pflag(b, 0, k, j, i)) || - (spherical && G.r(k, j, i) < r_eh + 0.1)) { // TODO customizable FOFC radius + (spherical && G.r(k, j, i) < r_eh + eh_buffer)) { fofcflag(b, 0, k, j, i) = 1; } else { fofcflag(b, 0, k, j, i) = 0; } } ); + + if (spherical && polar_cells > 0) { + for (int i_block = 0; i_block < guess->NumBlocks(); i_block++) { + auto &rc = guess->GetBlockData(i_block); + auto pmb = rc->GetBlockPointer(); + const bool is_inner_x2 = pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user; + const bool is_outer_x2 = pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user; + if (is_inner_x2 || is_outer_x2) { + auto lfofcflag = rc->PackVariables(std::vector{"fofcflag"}); + if (is_inner_x2) { + const IndexRange3 b = KDomain::GetRange(guess, IndexDomain::inner_x2); + int jstart = b.je + 1; + int jend = jstart + polar_cells - 1; + pmb0->par_for("fofc_mark_inner_x2", b.ks, b.ke, jstart, jend, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + lfofcflag(0, k, j, i) = 1; + } + ); + } + if (is_outer_x2) { + const IndexRange3 b = KDomain::GetRange(guess, IndexDomain::outer_x2); + int jend = b.js - 1; + int jstart = jend - polar_cells + 1; + pmb0->par_for("fofc_mark_outer_x2", b.ks, b.ke, jstart, jend, b.is, b.ie, + KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { + lfofcflag(0, k, j, i) = 1; + } + ); + } + } + } + } + return TaskStatus::complete; } @@ -108,16 +144,16 @@ TaskStatus Flux::FOFC(MeshData *md, MeshData *guess) const auto& U_all = md->PackVariablesAndFluxes(cons_flags, cons_map); const VarMap m_u(cons_map, true), m_p(prims_map, false); const int nvar = U_all.GetDim(4); + // Okay if this is empty since we won't access it then + const auto& Bf = md->PackVariables(std::vector{"cons.fB"}); // Parameters - const Real gam = pmb0->packages.Get("GRMHD")->Param("gamma"); - const bool use_global = pmb0->packages.Get("Flux")->Param("fofc_use_glf"); + const auto& pars = packages.Get("Flux")->AllParams(); + const Real gam = packages.Get("GRMHD")->Param("gamma"); + const bool use_global = pars.Get("fofc_use_glf"); const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(packages); - // With B_CT this will load the preference, without it will set face_b false and never access (empty) Bf - // Weird but it works - const bool face_b = (packages.AllPackages().count("B_CT") && - packages.Get("Flux")->Param("fofc_consistent_face_b")); - const auto& Bf = md->PackVariables(std::vector{"cons.fB"}); + // Only fix faces if they exist + const bool face_b = (Bf.GetDim(4) > 0 && pars.Get("fofc_consistent_face_b")); for (int dir=1; dir <= ndim; dir++) { // TODO if(trivial_direction) etc const TE el = FaceOf(dir); diff --git a/kharma/flux/get_flux.hpp b/kharma/flux/get_flux.hpp index fd02811e..6e482a75 100644 --- a/kharma/flux/get_flux.hpp +++ b/kharma/flux/get_flux.hpp @@ -79,13 +79,16 @@ inline TaskStatus GetFlux(MeshData *md) const bool reconstruction_floors = pars.Get("reconstruction_floors"); Floors::Prescription floors_temp; + Floors::Prescription floors_inner_temp; if (reconstruction_floors) { // Apply post-reconstruction floors. // Only enabled for WENO since it is not TVD, and only when other // floors are enabled. - floors_temp = packages.Get("Floors")->Param("prescription"); + floors_temp = packages.Get("Floors")->Param("prescription"); + floors_inner_temp = packages.Get("Floors")->Param("prescription_inner"); } const Floors::Prescription& floors = floors_temp; + const Floors::Prescription& floors_inner = floors_inner_temp; const bool reconstruction_fallback = pars.Get("reconstruction_fallback"); @@ -172,8 +175,8 @@ inline TaskStatus GetFlux(MeshData *md) // we have no guarantee they remotely resemble the *centered* primitives // If we selected to fall back to TVD, the floors are at zero (as intended) if (reconstruction_floors || reconstruction_fallback) { - fallback_tvd(i) = Floors::apply_geo_floors(G, Pl, m_p, gam, j, i, floors, loc); - fallback_tvd(i) |= Floors::apply_geo_floors(G, Pr, m_p, gam, j, i, floors, loc); + fallback_tvd(i) = Floors::apply_geo_floors(G, Pl, m_p, gam, j, i, floors, floors_inner, loc); + fallback_tvd(i) |= Floors::apply_geo_floors(G, Pr, m_p, gam, j, i, floors, floors_inner, loc); } } ); diff --git a/kharma/grmhd/grmhd.cpp b/kharma/grmhd/grmhd.cpp index b9e094fa..eac11bf9 100644 --- a/kharma/grmhd/grmhd.cpp +++ b/kharma/grmhd/grmhd.cpp @@ -94,7 +94,7 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

GetOrAddBoolean("parthenon/time", "start_dt_light", false); params.Add("start_dt_light", start_dt_light); bool use_dt_light = pin->GetOrAddBoolean("parthenon/time", "use_dt_light", false); @@ -110,6 +110,9 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

("type") == DriverType::imex) && (pin->GetBoolean("emhd", "on") || pin->GetOrAddBoolean("GRMHD", "implicit", false)); params.Add("implicit", implicit_grmhd); + // Explicitly-evolved ideal MHD variables as guess for Extended MHD runs + const bool ideal_guess = pin->GetOrAddBoolean("emhd", "ideal_guess", false); + params.Add("ideal_guess", ideal_guess); // AMR PARAMETERS // Adaptive mesh refinement options @@ -141,6 +144,12 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

>("cons_flags"); flags_cons.insert(flags_cons.end(), flags_grmhd.begin(), flags_grmhd.end()); + // Mark whether the ideal MHD variables are to be updated explicitly for the guess to the solver + if (ideal_guess) { + flags_prim.push_back(Metadata::GetUserFlag("IdealGuess")); + flags_cons.push_back(Metadata::GetUserFlag("IdealGuess")); + } + // We must additionally save the primtive variables as the "seed" for the next U->P solve flags_prim.push_back(Metadata::Restart); @@ -419,4 +428,55 @@ TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData *md) return TaskStatus::complete; } +void CancelBoundaryU3(MeshBlockData *rc, IndexDomain domain, bool coarse) +{ + // We're sometimes called on coarse buffers with or without AMR. + // Use of transmitting polar conditions when coarse buffers matter (e.g., refinement + // boundary touching the pole) is UNSUPPORTED + if (coarse) return; + + // Pull boundary properties + auto pmb = rc->GetBlockPointer(); + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const bool binner = KBoundaries::BoundaryIsInner(bface); + const auto bname = KBoundaries::BoundaryName(bface); + + // Pull variables (TODO take packs & maps, see boundaries.cpp) + PackIndexMap prims_map, cons_map; + auto P = rc->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map); + auto U = rc->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; + + const Real gam = pmb->packages.Get("GRMHD")->Param("gamma"); + + // Subtract the average B3 as "reconnection" + IndexRange3 b = KDomain::GetRange(rc, domain, coarse); + IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, coarse); + const int jf = (binner) ? bi.js : bi.je; // j index of last zone next to pole + parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_U3_" + bname, pmb->exec_space, + 0, 1, b.is, b.ie, + KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int& i) { + // Sum the first rank of U3 + Real U3_sum = 0.; + Kokkos::Sum sum_reducer(U3_sum); + parthenon::par_reduce_inner(member, bi.ks, bi.ke, + [&](const int& k, Real& local_result) { + local_result += isnan(P(m_p.U3, k, jf, i)) ? 0. : P(m_p.U3, k, jf, i); + } + , sum_reducer); + + // Calculate the average and subtract a portion + const Real U3_avg = U3_sum / (bi.ke - bi.ks + 1); + parthenon::par_for_inner(member, b.ks, b.ke, + [&](const int& k) { + P(m_p.U3, k, jf, i) -= U3_avg; + p_to_u(G, P, m_p, gam, k, jf, i, U, m_u); + } + ); + } + ); +} + } // namespace GRMHD diff --git a/kharma/grmhd/grmhd.hpp b/kharma/grmhd/grmhd.hpp index d9be5345..c528196a 100644 --- a/kharma/grmhd/grmhd.hpp +++ b/kharma/grmhd/grmhd.hpp @@ -88,4 +88,11 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin); * Currently just looks for negative density/internal energy */ TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData *rc); + +/** + * Subtract the average phi-component of velocity (U3) next to the pole, + * simulating polar velocity advecting together + */ +void CancelBoundaryU3(MeshBlockData *rc, IndexDomain domain, bool coarse); + } diff --git a/kharma/inverter/fixup.cpp b/kharma/inverter/fixup.cpp index 465cefae..c332475f 100644 --- a/kharma/inverter/fixup.cpp +++ b/kharma/inverter/fixup.cpp @@ -121,6 +121,9 @@ TaskStatus Inverter::FixUtoP(MeshBlockData *rc) const Floors::Prescription floors = pmb->packages.AllPackages().count("Floors") ? pmb->packages.Get("Floors")->Param("prescription") : pmb->packages.Get("Inverter")->Param("inverter_prescription"); + const Floors::Prescription floors_inner = pmb->packages.AllPackages().count("Floors") ? + pmb->packages.Get("Floors")->Param("prescription_inner") : + pmb->packages.Get("Inverter")->Param("inverter_prescription"); // We need the full packs of prims/cons for p_to_u // Pack new variables @@ -139,7 +142,7 @@ TaskStatus Inverter::FixUtoP(MeshBlockData *rc) if (failed(pflag(k, j, i))) { // Make sure all fixed values still abide by floors // TODO Full floors instead of just geo? - Floors::apply_geo_floors(G, P, m_p, gam, k, j, i, floors); + Floors::apply_geo_floors(G, P, m_p, gam, k, j, i, floors, floors_inner); // Make sure to keep lockstep // This will only be run for GRMHD, so we can call its p_to_u diff --git a/kharma/inverter/invert_template.hpp b/kharma/inverter/invert_template.hpp index ec3f7268..5026e728 100644 --- a/kharma/inverter/invert_template.hpp +++ b/kharma/inverter/invert_template.hpp @@ -69,6 +69,13 @@ KOKKOS_INLINE_FUNCTION bool failed(T status_flag) return static_cast(status_flag) > static_cast(Status::success); } +template +KOKKOS_INLINE_FUNCTION bool valid(T status_flag) +{ + // Return only values >0, among the failure flags + return static_cast(status_flag) == static_cast(Status::success); +} + /** * Recover local primitive variables, with a one-dimensional Newton-Raphson iterative solver. * Iteration starts from the current primitive values, and otherwse may *fail to converge* @@ -87,6 +94,6 @@ template KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const VariablePack& U, const VarMap& m_u, const Real& gam, const int& k, const int& j, const int& i, const VariablePack& P, const VarMap& m_p, - const Loci& loc, const Floors::Prescription& inverter_floors, + const Loci& loc, const Floors::Prescription& floors, const int& max_iterations, const Real& tol); } // namespace Inverter diff --git a/kharma/inverter/inverter.cpp b/kharma/inverter/inverter.cpp index db11afb7..7edc0008 100644 --- a/kharma/inverter/inverter.cpp +++ b/kharma/inverter/inverter.cpp @@ -68,8 +68,13 @@ std::shared_ptr Inverter::Initialize(ParameterInput *pin, std::sh // Use a custom block for inverter floors to allow customization. Not sure anyone *wants* that but... if (!pin->DoesBlockExist("inverter_floors")) { params.Add("inverter_prescription", Floors::MakePrescription(pin, "floors")); + if (pin->DoesBlockExist("floors_inner")) + params.Add("inverter_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "floors"), "floors_inner")); + else + params.Add("inverter_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "floors"), "floors")); } else { params.Add("inverter_prescription", Floors::MakePrescription(pin, "inverter_floors")); + params.Add("inverter_prescription_inner", Floors::MakePrescriptionInner(pin, Floors::MakePrescription(pin, "inverter_floors"), "inverter_floors")); } // Fixup options @@ -114,7 +119,6 @@ template inline void BlockPerformInversion(MeshBlockData *rc, IndexDomain domain, bool coarse) { auto pmb = rc->GetBlockPointer(); - const auto& G = pmb->coords; PackIndexMap prims_map, cons_map; auto U = GRMHD::PackMHDCons(rc, cons_map); @@ -132,18 +136,25 @@ inline void BlockPerformInversion(MeshBlockData *rc, IndexDomain domain, b auto &pars = pmb->packages.Get("Inverter")->AllParams(); const Real err_tol = pars.Get("err_tol"); const int iter_max = pars.Get("iter_max"); - Floors::Prescription inverter_floors = pars.Get("inverter_prescription"); + const Floors::Prescription inverter_floors = pars.Get("inverter_prescription"); + const Floors::Prescription inverter_floors_inner = pars.Get("inverter_prescription_inner"); + const bool radius_dependent_floors = inverter_floors.radius_dependent_floors; + + const auto& G = pmb->coords; // Get the primitives from our conserved versions // Notice 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); - 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 + && G.coords.is_spherical() + && G.r(k, j, i) < inverter_floors.floors_switch_r) ? + inverter_floors_inner : inverter_floors; int pflagl = Inverter::u_to_p(G, U, m_u, gam, k, j, i, P, m_p, Loci::center, - inverter_floors, iter_max, err_tol); + myfloors, iter_max, err_tol); pflag(0, k, j, i) = pflagl % Floors::FFlag::MINIMUM; int fflagl = (pflagl / Floors::FFlag::MINIMUM) * Floors::FFlag::MINIMUM; fflag(0, k, j, i) = fflagl; diff --git a/kharma/inverter/kastaun.hpp b/kharma/inverter/kastaun.hpp index 7834527e..3e536a21 100644 --- a/kharma/inverter/kastaun.hpp +++ b/kharma/inverter/kastaun.hpp @@ -215,7 +215,7 @@ template <> KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const VariablePack& U, const VarMap& m_u, const Real& gam, const int& k, const int& j, const int& i, const VariablePack& P, const VarMap& m_p, - const Loci& loc, const Floors::Prescription& inverter_floors, + const Loci& loc, const Floors::Prescription& floors, const int& max_iterations, const Real& tol) { // Shouldn't need this, KHARMA should die on NaN @@ -228,7 +228,8 @@ KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const V const Real a_over_g = alpha / G.gdet(loc, j, i); const Real D = U(m_u.RHO, k, j, i) * a_over_g; - const Real D_fl = std::max(D, inverter_floors.rho_min_const); + + const Real D_fl = std::max(D, floors.rho_min_const); Real Qcov[GR_DIM] = {(U(m_u.UU, k, j, i) - U(m_u.RHO, k, j, i)) * a_over_g, U(m_u.U1, k, j, i) * a_over_g, @@ -289,13 +290,12 @@ KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const V } //const Real zsq = rsq / h0sq_; // h0sq_ normalization set to 1 in Phoebus const Real zsq = rsq; - const Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / SQR(inverter_floors.gamma_max)); + const Real v0sq = std::min(zsq / (1.0 + zsq), 1.0 - 1.0 / SQR(floors.gamma_max)); // residual object. Caches most arguments/floors so calls are single-argument KastaunResidual res(D, q, bsq, bsq_rpsq, rsq, rbsq, v0sq, gam, - inverter_floors.rho_min_const, - inverter_floors.u_min_const/D_fl, - inverter_floors.gamma_max, inverter_floors.u_over_rho_max); + floors.rho_min_const, floors.u_min_const / D_fl, + floors.gamma_max, floors.u_over_rho_max); // SOLVE // TODO(BSP) better or faster solver? (Optionally) skip bracketing? diff --git a/kharma/inverter/onedw.hpp b/kharma/inverter/onedw.hpp index e3402d92..42c2d233 100644 --- a/kharma/inverter/onedw.hpp +++ b/kharma/inverter/onedw.hpp @@ -90,11 +90,12 @@ template <> KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const VariablePack& U, const VarMap& m_u, const Real& gam, const int& k, const int& j, const int& i, const VariablePack& P, const VarMap& m_p, - const Loci& loc, const Floors::Prescription& inverter_floors, + const Loci& loc, const Floors::Prescription& floors, const int& max_iterations, const Real& tol) { // TODO try inline floors in the old 1Dw? Probably not relevant anymore // Catch negative density + // TODO(BSP) if we start to see this hit, fix it here by returning P=P_floors if (U(m_u.RHO, k, j, i) <= 0.) { return static_cast(Status::neg_input); } @@ -130,6 +131,11 @@ KOKKOS_INLINE_FUNCTION int u_to_p(const GRCoordinates& G, const Var const Real QdB = dot(Bcon, Qcov); const Real Qdotn = dot(Qcon, ncov); + // TODO(BSP) check, test if this gets hit unlike above + // if (U(m_u.UU, k, j, i) <= gdet*(1e-8) + gdet*Bsq/2) { + + // } + Real Qtcon[GR_DIM]; DLOOP1 Qtcon[mu] = Qcon[mu] + ncon[mu] * Qdotn; const Real Qtsq = dot(Qcon, Qcov) + Qdotn*Qdotn; diff --git a/kharma/kharma.cpp b/kharma/kharma.cpp index 0e257eed..618a7de6 100644 --- a/kharma/kharma.cpp +++ b/kharma/kharma.cpp @@ -261,6 +261,40 @@ void KHARMA::FixParameters(ParameterInput *pin) if (tmp_coords.stopx(3) >= 0) pin->GetOrAddReal("parthenon/mesh", "x3max", tmp_coords.stopx(3)); + // Also set x1 refinements as a proportion of size + // TODO all regions! + if (pin->DoesBlockExist("parthenon/static_refinement0")) { + Real startx1 = pin->GetReal("parthenon/mesh", "x1min"); + Real stopx1 = pin->GetReal("parthenon/mesh", "x1max"); + Real lx1 = stopx1 - startx1; + Real startx1_prop = pin->GetReal("parthenon/static_refinement0", "x1min"); + Real stopx1_prop = pin->GetReal("parthenon/static_refinement0", "x1max"); + //std::cerr << "StartX1 " << startx1 << " lx1 " << lx1 << "Prop " << startx1_prop << " " << stopx1_prop << std::endl; + //std::cerr << "Adjust X1 " << startx1_prop*lx1 + startx1 << " to " << stopx1_prop*lx1 + startx1 << std::endl; + pin->SetReal("parthenon/static_refinement0", "x1min", std::max(startx1_prop*lx1 + startx1, startx1)); + pin->SetReal("parthenon/static_refinement0", "x1max", std::min(stopx1_prop*lx1 + startx1, stopx1)); + + if (pin->DoesParameterExist("parthenon/static_refinement0", "x2min")) { + Real startx2 = pin->GetReal("parthenon/mesh", "x2min"); + Real stopx2 = pin->GetReal("parthenon/mesh", "x2max"); + Real lx2 = stopx2 - startx2; + Real startx2_prop = pin->GetReal("parthenon/static_refinement0", "x2min"); + Real stopx2_prop = pin->GetReal("parthenon/static_refinement0", "x2max"); + pin->SetReal("parthenon/static_refinement0", "x2min", std::max(startx2_prop*lx2 + startx2, startx2)); + pin->SetReal("parthenon/static_refinement0", "x2max", std::min(stopx2_prop*lx2 + startx2, stopx2)); + } + + if (pin->DoesParameterExist("parthenon/static_refinement0", "x3min")) { + Real startx3 = pin->GetReal("parthenon/mesh", "x3min"); + Real stopx3 = pin->GetReal("parthenon/mesh", "x3max"); + Real lx3 = stopx3 - startx3; + Real startx3_prop = pin->GetReal("parthenon/static_refinement0", "x3min"); + Real stopx3_prop = pin->GetReal("parthenon/static_refinement0", "x3max"); + pin->SetReal("parthenon/static_refinement0", "x3min", std::max(startx3_prop*lx3 + startx3, startx3)); + pin->SetReal("parthenon/static_refinement0", "x3max", std::min(stopx3_prop*lx3 + startx3, stopx3)); + } + } + EndFlag(); } @@ -300,15 +334,22 @@ Packages_t KHARMA::ProcessPackages(std::unique_ptr &pin) // GRMHD needs globals to mark packages auto t_grmhd = tl.AddTask(t_globals | t_driver, KHARMA::AddPackage, packages, GRMHD::Initialize, pin.get()); // Only load the inverter if GRMHD/EMHD isn't being evolved implicitly + // Unless we want to use the explicitly-evolved ideal MHD variables as a guess for the solver auto t_inverter = t_grmhd; - if (!pin->GetOrAddBoolean("GRMHD", "implicit", pin->GetOrAddBoolean("emhd", "on", false))) { + if (!pin->GetOrAddBoolean("GRMHD", "implicit", pin->GetOrAddBoolean("emhd", "on", false)) || + pin->GetOrAddBoolean("emhd", "ideal_guess", false)) { t_inverter = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, Inverter::Initialize, pin.get()); } - // Floors package is only loaded if floors aren't disabled (TODO rename "on"?) - if (!pin->GetOrAddBoolean("floors", "disable_floors", false)) { + // Floors package is only loaded if floors aren't disabled + // Respect legacy version for a while + bool floors_on_default = true; + if (pin->DoesParameterExist("floors", "disable_floors")) { + floors_on_default = !pin->GetBoolean("floors", "disable_floors"); + } + if (pin->GetOrAddBoolean("floors", "on", floors_on_default)) { auto t_floors = tl.AddTask(t_inverter, KHARMA::AddPackage, packages, Floors::Initialize, pin.get()); } - // Reductions, needed for most other packages + // Reductions, needed by most other packages auto t_reductions = tl.AddTask(t_none, KHARMA::AddPackage, packages, Reductions::Initialize, pin.get()); // B field solvers, to ensure divB ~= 0. @@ -324,7 +365,7 @@ Packages_t KHARMA::ProcessPackages(std::unique_ptr &pin) } else if (b_field_solver == "constrained_transport" || b_field_solver == "face_ct") { t_b_field = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, B_CT::Initialize, pin.get()); } else if (b_field_solver == "constraint_damping" || b_field_solver == "cd") { - // Constraint damping, probably only useful for non-GR MHD systems + // Constraint damping. NON-WORKING t_b_field = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, B_CD::Initialize, pin.get()); } else if (b_field_solver == "flux_ct") { t_b_field = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, B_FluxCT::Initialize, pin.get()); @@ -353,7 +394,7 @@ Packages_t KHARMA::ProcessPackages(std::unique_ptr &pin) if (pin->GetOrAddBoolean("electrons", "on", false)) { auto t_electrons = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, Electrons::Initialize, pin.get()); } - if (pin->GetBoolean("emhd", "on")) { + if (pin->GetBoolean("emhd", "on")) { // Set above when deciding to load inverter auto t_emhd = tl.AddTask(t_grmhd, KHARMA::AddPackage, packages, EMHD::Initialize, pin.get()); } if (pin->GetOrAddBoolean("wind", "on", false)) { diff --git a/kharma/kharma_package.hpp b/kharma/kharma_package.hpp index 3d53e181..875d7b82 100644 --- a/kharma/kharma_package.hpp +++ b/kharma/kharma_package.hpp @@ -159,8 +159,7 @@ TaskStatus AddSource(MeshData *md, MeshData *mdudt, IndexDomain doma TaskStatus MeshApplyPrimSource(MeshData *md); /** - * Apply all floors, including any package-specific limiters. - * This function respects "disable_floors". + * Apply all registered floors, including any package-specific limiters. * * LOCKSTEP: this function respects P and returns consistent P<->U */ diff --git a/kharma/prob/bondi.cpp b/kharma/prob/bondi.cpp index 4165e77f..4d159d34 100644 --- a/kharma/prob/bondi.cpp +++ b/kharma/prob/bondi.cpp @@ -38,13 +38,8 @@ #include "floors.hpp" #include "flux_functions.hpp" -/** - * Initialization of a Bondi problem with specified sonic point & accretion rate - */ -TaskStatus InitializeBondi(std::shared_ptr>& rc, ParameterInput *pin) +void AddBondiParameters(ParameterInput *pin, Packages_t &packages) { - auto pmb = rc->GetBlockPointer(); - const Real mdot = pin->GetOrAddReal("bondi", "mdot", 1.0); const Real rs = pin->GetOrAddReal("bondi", "rs", 8.0); const Real ur_frac = pin->GetOrAddReal("bondi", "ur_frac", 1.); @@ -69,22 +64,37 @@ TaskStatus InitializeBondi(std::shared_ptr>& rc, ParameterIn // Add these to package properties, since they continue to be needed on boundaries // TODO Problems NEED params - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("mdot")) - pmb->packages.Get("GRMHD")->AddParam("mdot", mdot); - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("rs")) - pmb->packages.Get("GRMHD")->AddParam("rs", rs); - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("rin_bondi")) - pmb->packages.Get("GRMHD")->AddParam("rin_bondi", rin_bondi); - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("fill_interior_bondi")) - pmb->packages.Get("GRMHD")->AddParam("fill_interior_bondi", fill_interior); - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("zero_velocity_bondi")) - pmb->packages.Get("GRMHD")->AddParam("zero_velocity_bondi", zero_velocity); - if(! pmb->packages.Get("GRMHD")->AllParams().hasKey("diffinit_bondi")) - pmb->packages.Get("GRMHD")->AddParam("diffinit_bondi", diffinit); - if(! (pmb->packages.Get("GRMHD")->AllParams().hasKey("ur_frac"))) - pmb->packages.Get("GRMHD")->AddParam("ur_frac", ur_frac); - if(! (pmb->packages.Get("GRMHD")->AllParams().hasKey("uphi"))) - pmb->packages.Get("GRMHD")->AddParam("uphi", uphi); + if(! packages.Get("GRMHD")->AllParams().hasKey("mdot")) + packages.Get("GRMHD")->AddParam("mdot", mdot); + if(! packages.Get("GRMHD")->AllParams().hasKey("rs")) + packages.Get("GRMHD")->AddParam("rs", rs); + if(! packages.Get("GRMHD")->AllParams().hasKey("rin_bondi")) + packages.Get("GRMHD")->AddParam("rin_bondi", rin_bondi); + if(! packages.Get("GRMHD")->AllParams().hasKey("fill_interior_bondi")) + packages.Get("GRMHD")->AddParam("fill_interior_bondi", fill_interior); + if(! packages.Get("GRMHD")->AllParams().hasKey("zero_velocity_bondi")) + packages.Get("GRMHD")->AddParam("zero_velocity_bondi", zero_velocity); + if(! packages.Get("GRMHD")->AllParams().hasKey("diffinit_bondi")) + packages.Get("GRMHD")->AddParam("diffinit_bondi", diffinit); + if(! (packages.Get("GRMHD")->AllParams().hasKey("ur_frac"))) + packages.Get("GRMHD")->AddParam("ur_frac", ur_frac); + if(! (packages.Get("GRMHD")->AllParams().hasKey("uphi"))) + packages.Get("GRMHD")->AddParam("uphi", uphi); +} + +/** + * Initialization of a Bondi problem with specified sonic point & accretion rate + */ +TaskStatus InitializeBondi(std::shared_ptr>& rc, ParameterInput *pin) +{ + auto pmb = rc->GetBlockPointer(); + + // Add parameters we'll need throughout the run to 'GRMHD' package + AddBondiParameters(pin, pmb->packages); + + // We need these two + const Real rin_bondi = pmb->packages.Get("GRMHD")->Param("rin_bondi"); + const bool fill_interior = pmb->packages.Get("GRMHD")->Param("fill_interior_bondi"); // Set this problem to control the outer X1 boundary by default // remember to disable inflow_check in parameter file! @@ -116,7 +126,7 @@ TaskStatus InitializeBondi(std::shared_ptr>& rc, ParameterIn // Apply floors to initialize the any part of the domain we didn't // Bondi's BL coordinates do not like the EH, so we replace the zeros with something reasonable - // Note this ignores the "disable_floors" parameter, since it's necessary for initialization + // Note this ignores whether the "Floors" package is loaded, since it's necessary for initialization if (rin_bondi > pin->GetReal("coordinates", "r_in") && !(fill_interior)) { Floors::ApplyInitialFloors(pin, rc.get(), IndexDomain::interior); } diff --git a/kharma/prob/bondi.hpp b/kharma/prob/bondi.hpp index 4ca0d4a4..43d3a05c 100644 --- a/kharma/prob/bondi.hpp +++ b/kharma/prob/bondi.hpp @@ -50,6 +50,12 @@ */ TaskStatus InitializeBondi(std::shared_ptr>& rc, ParameterInput *pin); +/** + * Record parameters Bondi problem (or boundaries!) will need throughout the run + * Currently uses "GRMHD" package as convenient proxy (TODO fix that with problem-packages) + */ +void AddBondiParameters(ParameterInput *pin, Packages_t &packages); + /** * Set all values on a given domain to the Bondi inflow analytic steady-state solution. * Use the template version when possible, which just calls through diff --git a/kharma/prob/bz_monopole.cpp b/kharma/prob/bz_monopole.cpp deleted file mode 100644 index 3a442091..00000000 --- a/kharma/prob/bz_monopole.cpp +++ /dev/null @@ -1,84 +0,0 @@ -/* - * File: bz_monopole.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 - * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR - * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER - * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, - * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -#include "bz_monopole.hpp" - -#include "coordinate_utils.hpp" -#include "types.hpp" - -#include -#include "Kokkos_Random.hpp" - -TaskStatus InitializeBZMonopole(std::shared_ptr>& rc, ParameterInput *pin) -{ - auto pmb = rc->GetBlockPointer(); - GridScalar rho = rc->Get("prims.rho").data; - GridScalar u = rc->Get("prims.u").data; - GridVector uvec = rc->Get("prims.uvec").data; - - Real bsq_o_rho_max = pin->GetOrAddReal("floors", "bsq_over_rho_max", 1.e2); - Real rho_min_limit = pin->GetOrAddReal("floors", "rho_min_geom", 1.e-6); - Real u_min_limit = pin->GetOrAddReal("floors", "u_min_geom", 1.e-8); - - 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); - int ks = pmb->cellbounds.ks(domain), ke = pmb->cellbounds.ke(domain); - - const auto& G = pmb->coords; - const GReal a = G.coords.get_a(); - - pmb->par_for("fm_torus_init", ks, ke, js, je, is, ie, - KOKKOS_LAMBDA (const int &k, const int &j, const int &i) { - GReal Xembed[GR_DIM]; - G.coord_embed(k, j, i, Loci::center, Xembed); - GReal r = Xembed[1]; - - GReal r_horizon = 1. + m::sqrt(1. - a*a); - GReal r_char = 10. * r_horizon; - - Real trho = rho_min_limit + (r / r_char) / m::pow(r, 4.) / bsq_o_rho_max; - Real tu = u_min_limit + (r / r_char) / m::pow(r, 4.) / bsq_o_rho_max; - - rho(k, j, i) = trho; - u(k, j, i) = tu; - uvec(0, k, j, i) = 0.; - uvec(1, k, j, i) = 0.; - uvec(2, k, j, i) = 0.; - } - ); - - return TaskStatus::complete; -} - diff --git a/kharma/prob/bz_monopole.hpp b/kharma/prob/bz_monopole.hpp deleted file mode 100644 index a0b69f72..00000000 --- a/kharma/prob/bz_monopole.hpp +++ /dev/null @@ -1,11 +0,0 @@ -// Monopole problem initialization with floors -#pragma once - -#include "decs.hpp" -#include "types.hpp" - -/** - * Initialize a Blandford-Znajek monopole setup - */ -TaskStatus InitializeBZMonopole(std::shared_ptr>& rc, ParameterInput *pin); - diff --git a/kharma/prob/fm_torus.cpp b/kharma/prob/fm_torus.cpp index 879c1b34..736fe68e 100644 --- a/kharma/prob/fm_torus.cpp +++ b/kharma/prob/fm_torus.cpp @@ -187,8 +187,8 @@ TaskStatus InitializeFMTorus(std::shared_ptr>& rc, Parameter } ); - // Apply floors to initialize the rest of the domain (regardless of the 'disable_floors' param) - // Since the conserved vars U are not initialized, this is done in *fluid frame*, + // Apply floors to initialize the rest of the domain (regardless whether we'll use them later) + // Since the conserved vars U are not initialized, this is effectively done in *fluid frame*, // even if NOF frame is chosen (iharm3d does the same iiuc) // This is probably not a huge issue, just good to state explicitly Floors::ApplyInitialFloors(pin, rc.get(), IndexDomain::interior); diff --git a/kharma/prob/gizmo.cpp b/kharma/prob/gizmo.cpp index 0e604f50..1572aff9 100644 --- a/kharma/prob/gizmo.cpp +++ b/kharma/prob/gizmo.cpp @@ -146,7 +146,7 @@ TaskStatus SetGIZMO(std::shared_ptr>& rc, IndexDomain domain int i_sh; GReal del_sh; XtoindexGIZMO(Xshell, r_device, length, i_sh, del_sh); - Real vacuum_rho, vacuum_u_over_rho, vacuum_logrho, vacuum_log_u_over_rho; + Real vacuum_rho, vacuum_u_over_rho; vacuum_rho = rho_device(i_sh)*(1.-del_sh)+rho_device(i_sh+1)*del_sh; vacuum_u_over_rho = (T_device(i_sh)*(1.-del_sh)+T_device(i_sh+1)*del_sh)/(gam-1.); diff --git a/kharma/prob/problem.cpp b/kharma/prob/problem.cpp index e00c0d69..cd9a85ca 100644 --- a/kharma/prob/problem.cpp +++ b/kharma/prob/problem.cpp @@ -51,7 +51,6 @@ #include "resize_restart.hpp" #include "resize_restart_kharma.hpp" #include "kelvin_helmholtz.hpp" -#include "bz_monopole.hpp" #include "mhdmodes.hpp" #include "orszag_tang.hpp" #include "shock_tube.hpp" diff --git a/kharma/prob/seed_B.cpp b/kharma/prob/seed_B.cpp index 0ef7b0f3..daff2cbe 100644 --- a/kharma/prob/seed_B.cpp +++ b/kharma/prob/seed_B.cpp @@ -112,7 +112,6 @@ TaskStatus SeedBFieldType(MeshBlockData *rc, ParameterInput *pin, IndexDom // Shortcut to field values for easy fields if constexpr (Seed == BSeedType::constant || Seed == BSeedType::monopole || - Seed == BSeedType::monopole_cube || Seed == BSeedType::orszag_tang || Seed == BSeedType::wave || Seed == BSeedType::shock_tube) @@ -425,8 +424,8 @@ TaskStatus SeedBField(MeshData *md, ParameterInput *pin) status = SeedBFieldType(rc, pin); } else if (b_field_type == "monopole") { status = SeedBFieldType(rc, pin); - } else if (b_field_type == "monopole_cube") { - status = SeedBFieldType(rc, pin); + } else if (b_field_type == "monopole_cube") { // Legacy name for the correct monopole init + status = SeedBFieldType(rc, pin); } else if (b_field_type == "sane") { status = SeedBFieldType(rc, pin); } else if (b_field_type == "mad") { diff --git a/kharma/prob/seed_B.hpp b/kharma/prob/seed_B.hpp index 71b1c1b6..41d6fa97 100644 --- a/kharma/prob/seed_B.hpp +++ b/kharma/prob/seed_B.hpp @@ -51,7 +51,7 @@ TaskStatus NormalizeBField(MeshData *md, ParameterInput *pin); */ // Internal representation of the field initialization preference, used for templating -enum BSeedType{constant, monopole, monopole_cube, orszag_tang, orszag_tang_a, wave, shock_tube, +enum BSeedType{constant, monopole, orszag_tang, orszag_tang_a, wave, shock_tube, sane, mad, mad_quadrupole, r3s3, r5s5, gaussian, bz_monopole, vertical, r1s2}; #define SEEDA_ARGS GReal *x, const GReal *dxc, double rho, double rin, double min_A, double A0, double arg1, double rb @@ -150,16 +150,9 @@ KOKKOS_INLINE_FUNCTION void seed_b(SEEDB_ARGS) { B1 = 0./0.; B2 = 0./0.; B3 = 0. template<> KOKKOS_INLINE_FUNCTION void seed_b(SEEDB_ARGS) {} -// Reduce radial component by gdet for constant flux -template<> -KOKKOS_INLINE_FUNCTION void seed_b(SEEDB_ARGS) -{ - B1 /= gdet; -} - // Reduce radial component by the cube of radius template<> -KOKKOS_INLINE_FUNCTION void seed_b(SEEDB_ARGS) +KOKKOS_INLINE_FUNCTION void seed_b(SEEDB_ARGS) { B1 /= (x[1]*x[1]*x[1]); } diff --git a/kharma/prob/shock_tube.hpp b/kharma/prob/shock_tube.hpp index b9f5abcf..b191bdaa 100644 --- a/kharma/prob/shock_tube.hpp +++ b/kharma/prob/shock_tube.hpp @@ -39,7 +39,7 @@ TaskStatus InitializeShockTube(std::shared_ptr>& rc, Paramet const Real B3L = pin->GetOrAddReal("shock", "B3L", 0.0); const Real B3R = pin->GetOrAddReal("shock", "B3R", 0.0); - IndexDomain domain = IndexDomain::interior; + IndexDomain domain = IndexDomain::entire; IndexRange ib = pmb->cellbounds.GetBoundsI(domain); IndexRange jb = pmb->cellbounds.GetBoundsJ(domain); IndexRange kb = pmb->cellbounds.GetBoundsK(domain); diff --git a/kharma/types.hpp b/kharma/types.hpp index 15d1678e..e35e2ef5 100644 --- a/kharma/types.hpp +++ b/kharma/types.hpp @@ -71,14 +71,20 @@ constexpr TE E3 = TE::E3; constexpr TE NN = TE::NN; // Any basic type manips, see LocOf in decs etc etc +// Host-only because the templating/types are weird on device template -KOKKOS_FORCEINLINE_FUNCTION TopologicalElement FaceOf(const T& dir) { +inline TE FaceOf(const T& dir) { return (dir == X1DIR) ? F1 : ((dir == X2DIR) ? F2 : F3); } template -KOKKOS_FORCEINLINE_FUNCTION TopologicalElement EdgeOf(const T& dir) { +inline TE EdgeOf(const T& dir) { return (dir == X1DIR) ? E1 : ((dir == X2DIR) ? E2 : E3); } +template +inline std::vector OrthogonalEdges(const T& dir) { + return (dir == X1DIR) ? std::vector{E2, E3} : + ((dir == X2DIR) ? std::vector{E1, E3} : std::vector{E1, E2}); +} // Struct for derived 4-vectors at a point, usually calculated and needed together typedef struct { @@ -89,12 +95,12 @@ typedef struct { } FourVectors; typedef struct { - uint is; - uint ie; - uint js; - uint je; - uint ks; - uint ke; + int is; + int ie; + int js; + int je; + int ks; + int ke; } IndexRange3; typedef struct { diff --git a/machines/darwin.sh b/machines/darwin.sh index 2d9b86cb..1591a1d2 100644 --- a/machines/darwin.sh +++ b/machines/darwin.sh @@ -5,7 +5,7 @@ if [[ ($HOSTNAME == "cn"* || $HOSTNAME == "darwin"*) && ("$PWD" == "/projects/jacamar-ci"* || "$PWD" == "/vast"*) ]]; then - module purge + #module purge module load cmake # Where we're going, we don't need system libraries @@ -32,9 +32,6 @@ if [[ ($HOSTNAME == "cn"* || $HOSTNAME == "darwin"*) && module load nvhpc C_NATIVE="nvc" CXX_NATIVE="nvc++" - # New NVHPC doesn't like CUDA_HOME - export NVHPC_CUDA_HOME="$CUDA_HOME" - unset CUDA_HOME elif [[ "$ARGS" == *"icc"* ]]; then module load intel-classic/2021.3.0 C_NATIVE=icc @@ -45,9 +42,6 @@ if [[ ($HOSTNAME == "cn"* || $HOSTNAME == "darwin"*) && module load nvhpc C_NATIVE="nvc" CXX_NATIVE="nvc++" - # New NVHPC doesn't like CUDA_HOME - export NVHPC_CUDA_HOME="$CUDA_HOME" - unset CUDA_HOME else module load intel C_NATIVE=icx @@ -57,8 +51,13 @@ if [[ ($HOSTNAME == "cn"* || $HOSTNAME == "darwin"*) && # 2. Load accelerator libraries if [[ "$ARGS" == *"cuda"* ]]; then - module load cuda/12.0.0 nvhpc - PREFIX_PATH=$NVHPC_ROOT + module load cuda/12.3.1 + # Newer NVHPC wants us to leave it alone + #unset CUDA_HOME + # For manually exporting CUDA and COMM_LIBS + #export NVHPC_CUDA_HOME="$CUDA_HOME" + #export NVHPC_COMM_LIBS_HOME=/projects/darwin-nv/rhel8/aarch64/packages/nvhpc/Linux_aarch64/24.1/comm_libs + #PREFIX_PATH=$NVHPC_ROOT #EXTRA_FLAGS="-DPARTHENON_ENABLE_HOST_COMM_BUFFERS=ON $EXTRA_FLAGS" elif [[ "$ARGS" == *"hip"* ]]; then # No MPI or OpenMP -- No OFI OpenMPI on Darwin (right?) and HIP hates OpenMP diff --git a/make.sh b/make.sh index d4b8a530..d81be542 100755 --- a/make.sh +++ b/make.sh @@ -174,10 +174,6 @@ elif [[ "$ARGS" == *"cuda"* ]]; then echo "Dry-running the nvcc wrapper with $CXXFLAGS" fi export NVCC_WRAPPER_DEFAULT_COMPILER="$CXX_NATIVE" - # Generally Kokkos sets this, so we don't need to - #export CXXFLAGS="--expt-relaxed-constexpr $CXXFLAGS" - # New NVHPC complains if we don't set this - export NVHPC_CUDA_HOME=$CUDA_HOME OUTER_LAYOUT="MANUAL1D_LOOP" INNER_LAYOUT="TVR_INNER_LOOP" ENABLE_OPENMP="ON" diff --git a/pars/smr/kelvin_helmholtz.par b/pars/smr/kelvin_helmholtz.par index 005657ec..6eb6ace6 100644 --- a/pars/smr/kelvin_helmholtz.par +++ b/pars/smr/kelvin_helmholtz.par @@ -35,10 +35,8 @@ nx3 = 1 x1min = 0.4 x1max = 0.6 -x2min = 0.9 -x2max = 1.1 -x3min = 0.0 -x3max = 0.0 +x2min = 0.45 +x2max = 0.55 level = 1 diff --git a/pars/smr/mhdmodes_refined.par b/pars/smr/mhdmodes_refined.par index 1b1901fa..4bf94a88 100644 --- a/pars/smr/mhdmodes_refined.par +++ b/pars/smr/mhdmodes_refined.par @@ -63,6 +63,7 @@ implicit = false solver = face_ct +ct_scheme = bs99 implicit = false diff --git a/pars/smr/orszag_tang_refined.par b/pars/smr/orszag_tang_refined.par index 4f41cf1f..9d7ac508 100644 --- a/pars/smr/orszag_tang_refined.par +++ b/pars/smr/orszag_tang_refined.par @@ -25,10 +25,14 @@ nx2 = 64 nx3 = 1 -x1min = -0.1 -x1max = 0.1 -x2min = -0.1 -x2max = 0.1 +# Refinement is expressed as an N-dimensional convex box, +# with each coordinate given as a proportion 0.0-1.0 of its dimension +# Any meshblock intersecting the box gets refined to the given level +# e.g. this is the central region ~0 +x1min = 0.49 +x1max = 0.51 +x2min = 0.49 +x2max = 0.51 level = 1 diff --git a/pars/smr/orszag_tang_refined_small.par b/pars/smr/orszag_tang_refined_small.par index 6f516862..0f9d4956 100644 --- a/pars/smr/orszag_tang_refined_small.par +++ b/pars/smr/orszag_tang_refined_small.par @@ -26,10 +26,14 @@ nx2 = 32 nx3 = 1 -x1min = -3.14 -x1max = -3.1 -x2min = -3.14 -x2max = -3.1 +# Refinement is expressed as an N-dimensional convex box, +# with each coordinate given as a proportion 0.0-1.0 of its dimension +# Any meshblock intersecting the box gets refined to the given level +# e.g. this is the bottom left corner, for testing +x1min = 0.0 +x1max = 0.1 +x2min = 0.0 +x2max = 0.1 level = 1 # Set boring box coordinates. Explanations in bondi.par diff --git a/pars/smr/sane2d_refined.par b/pars/smr/sane2d_refined.par index 729153aa..36c3b873 100644 --- a/pars/smr/sane2d_refined.par +++ b/pars/smr/sane2d_refined.par @@ -17,8 +17,15 @@ nx2 = 64 nx3 = 1 -x1min = 1.0 -x1max = 3.0 +# Refinement is expressed as an N-dimensional convex box, +# with each coordinate given as a proportion 0.0-1.0 of its dimension +# Any meshblock intersecting the box gets refined to the given level +# 5 base blocks +# + 3 blocks from refining the central one +# + 12 blocks from refining the 4 central blocks again +# = 20 blocks +x1min = 0.0 +x1max = 0.49 x2min = 0.49 x2max = 0.51 level = 2 diff --git a/pars/smr/sane3d_refined.par b/pars/smr/sane3d_refined.par index 2b3ee29e..3c634ca3 100644 --- a/pars/smr/sane3d_refined.par +++ b/pars/smr/sane3d_refined.par @@ -8,7 +8,7 @@ problem_id = torus refinement = static numlevel = 2 nx1 = 128 -nx2 = 64 +nx2 = 96 nx3 = 64 @@ -17,12 +17,17 @@ nx2 = 32 nx3 = 32 -x1min = 1.0 -x1max = 2.0 +# Refinement is expressed as an N-dimensional convex box, +# with each coordinate given as a proportion 0.0-1.0 of its dimension +# Any meshblock intersecting the box gets refined to the given level +# This refines the inner two radial belts of blocks -- +# 4*2*3 base + 2*1*2*7 additional from refinement = 52 blocks total +x1min = 0.0 +x1max = 0.49 x2min = 0.49 x2max = 0.51 x3min = 0.0 -x3max = 6.28 +x3max = 1.0 level = 1 diff --git a/scripts/ci/darwin.yml b/scripts/ci/darwin.yml index 9f7b4047..664567a5 100644 --- a/scripts/ci/darwin.yml +++ b/scripts/ci/darwin.yml @@ -57,7 +57,7 @@ build: - echo "Skipping pyharm install in build." script: - export PREFIX_PATH=$PWD/external/hdf5 - - ./make.sh clean cuda hdf5 volta ompi + - ./make.sh clean cuda hdf5 volta artifacts: paths: - kharma.* diff --git a/tests/bondi_viscous/run.sh b/tests/bondi_viscous/run.sh index a50a00bc..9759e022 100755 --- a/tests/bondi_viscous/run.sh +++ b/tests/bondi_viscous/run.sh @@ -22,7 +22,7 @@ conv_2d() { mv bondi.out0.final.phdf emhd_2d_${res}_end_${1}.phdf done check_code=0 - python check.py $ALL_RES $1 2d || check_code=$? + python3 check.py $ALL_RES $1 2d || check_code=$? rm -r *.xdmf rm -r *.out0* if [[ $check_code != 0 ]]; then @@ -34,6 +34,8 @@ conv_2d() { } ALL_RES="8,16,32,64" -conv_2d emhd2d_weno driver/reconstruction=weno5 "in 2D, WENO5" +conv_2d emhd2d_weno driver/reconstruction=weno5 "Viscous Bondi in 2D, WENO5" +# Test if it works with ideal solution as guess +conv_2d emhd2d_weno_ideal_guess emhd/ideal_guess=true "Viscous bondi in 2D, WENO5, Ideal guess" exit $exit_code diff --git a/tests/bz_monopole/bz_monopole.par b/tests/bz_monopole/bz_monopole.par index a0b26367..0f6071d8 100644 --- a/tests/bz_monopole/bz_monopole.par +++ b/tests/bz_monopole/bz_monopole.par @@ -18,13 +18,10 @@ nx3 = 1 base = spherical_ks -transform = fmks +transform = mks r_out = 100. a = 0.9375 hslope = 0.3 -mks_smooth = 0.5 -poly_xt = 0.82 -poly_alpha = 14.0 tlim = 100.0 @@ -33,22 +30,21 @@ nlim = -1 verbose = 1 extra_checks = 1 -flag_verbose = 0 +flag_verbose = 2 cfl = 0.7 gamma = 1.444444 -reconstruction = linear_mc +reconstruction = weno5 type = bz_monopole norm = false -bsq_over_rho_max = 100 -rho_min_geom = 1e-20 -u_min_geom = 1e-20 -gamma_max = 10 +rho_min_geom = 1e-6 +u_min_geom = 1e-8 +#bsq_over_rho_max = 1000 on = false @@ -61,7 +57,7 @@ power = 40 file_type = hdf5 dt = 5.0 single_precision_output = false -variables = prims.rho, prims.u, prims.uvec, prims.B, cons.B, divB +variables = prims, cons.B, divB ghost_zones = true diff --git a/tests/conducting_atmosphere/run.sh b/tests/conducting_atmosphere/run.sh index f1e801d9..8f18c3b7 100755 --- a/tests/conducting_atmosphere/run.sh +++ b/tests/conducting_atmosphere/run.sh @@ -25,7 +25,7 @@ conv_2d() { done check_code=0 pyharm-convert --double *.phdf - python check.py $ALL_RES $1 2d || check_code=$? + python3 check.py $ALL_RES $1 2d || check_code=$? if [[ $check_code != 0 ]]; then echo Conducting atmosphere test $3 FAIL: $check_code exit_code=1 @@ -36,5 +36,7 @@ conv_2d() { ALL_RES="64,128,256,512" conv_2d emhd2d_weno driver/reconstruction=weno5 "in 2D, WENO5" +# Test if it works with ideal solution as guess +conv_2d emhd2d_weno_ideal_guess emhd/ideal_guess=true "in 2D, WENO5, Ideal guess" exit $exit_code diff --git a/tests/emhdmodes/check.py b/tests/emhdmodes/check.py index 7c577498..472bcebc 100644 --- a/tests/emhdmodes/check.py +++ b/tests/emhdmodes/check.py @@ -97,10 +97,11 @@ if not (dvar_cos[k] == 0 and dvar_sin[k] == 0): powerfits[k] = np.polyfit(np.log(RES), np.log(L1[:,k]), 1)[0] print("Power fit {}: {} {}".format(VARS[k], powerfits[k], L1[:,k])) - if powerfits[k] > -1.6 or powerfits[k] < -2.7: + if powerfits[k] > -1.6 or powerfits[k] < -3.3: # Everything *should* converge at ~2, but we relax the reqt due to known behavior: # 1. B field in WENO seems to lag, at ~1.7 # 2. Problems run under linear/MC seem to converge ~2.5 in most variables + # 3. EMHD modes with ideal guess has ~3 convergence for rho fail = 1 # plot diff --git a/tests/emhdmodes/run.sh b/tests/emhdmodes/run.sh index b720c664..00b6e613 100755 --- a/tests/emhdmodes/run.sh +++ b/tests/emhdmodes/run.sh @@ -21,7 +21,7 @@ conv_2d() { mv emhdmodes.out0.final.phdf emhd_2d_${res}_end_${1}.phdf done check_code=0 - python check.py $ALL_RES "$3" $1 2d || check_code=$? + python3 check.py $ALL_RES "$3" $1 2d || check_code=$? if [[ $check_code != 0 ]]; then echo $3 FAIL: $check_code exit_code=1 @@ -34,12 +34,14 @@ conv_2d() { ALL_RES="32,64,128" conv_2d emhd2d_weno flux/reconstruction=weno5 "EMHD mode in 2D, WENO5" -# ...but linear doesn't capture wave until higher res. Troubling. +# ...but linear doesn't capture wave until higher res. Troubling. ALL_RES="64,128,256" conv_2d emhd2d_mc flux/reconstruction=linear_mc "EMHD mode in 2D, linear/MC reconstruction" # Test that higher-order terms don't mess anything up conv_2d emhd2d_higher_order emhd/higher_order_terms=true "EMHD mode in 2D, higher order terms enabled" # Test we can use imex/EMHD and face CT conv_2d emhd2d_face_ct b_field/solver=face_ct "EMHD mode in 2D w/Face CT" +# Test if it works with ideal solution as guess +conv_2d emhd2d_ideal_guess emhd/ideal_guess=true "EMHD mode in 2D, Ideal guess" exit $exit_code diff --git a/tests/mhdmodes/check.py b/tests/mhdmodes/check.py index 0dfffc21..9599cf50 100644 --- a/tests/mhdmodes/check.py +++ b/tests/mhdmodes/check.py @@ -141,16 +141,18 @@ fail = 1 # MAKE PLOTS -fig = plt.figure(figsize=(5,5)) +fig = plt.figure(figsize=(3,3)) ax = fig.add_subplot(1,1,1) for k in range(NVAR): if abs(dvar[k]) != 0.: ax.plot(RES, L1[:,k], marker='s', label=VARS[k]) -norm = L1[0,0]*RES[0]*RES[0] -if norm < 1e-4: - norm = L1[0,3]*RES[0]*RES[0] +if "alfven" in SHORT: + norm = L1[0,4]*RES[0]*RES[0] +else: + norm = L1[0,0]*RES[0]*RES[0] + xmin = RES[0]/2. xmax = RES[-1]*2. ax.plot([xmin, xmax], norm*np.asarray([xmin, xmax])**-2., color='k', linestyle='--', label='N^-2') @@ -158,8 +160,9 @@ plt.xscale('log', base=2); plt.yscale('log') plt.xlim([RES[0]/np.sqrt(2.), RES[-1]*np.sqrt(2.)]) plt.xlabel('N'); plt.ylabel('L1') -plt.title("{}".format(LONG)) +#plt.title("{}".format(LONG)) plt.legend(loc=1) -plt.savefig("convergence_modes_{}_{}.png".format(DIM,SHORT)) +plt.subplots_adjust(left=0.25, bottom=0.18, top=0.97, right=0.97) +plt.savefig("convergence_modes_{}_{}.png".format(DIM,SHORT), dpi=300) exit(fail) diff --git a/tests/mhdmodes/run.sh b/tests/mhdmodes/run.sh index 34802fd6..9b2fbdf5 100755 --- a/tests/mhdmodes/run.sh +++ b/tests/mhdmodes/run.sh @@ -29,7 +29,7 @@ conv_3d() { mv mhdmodes.out0.final.phdf mhd_3d_${1}_${res}_end.phdf done check_code=0 - python check.py $ALL_RES "$3" $1 3d 3 || check_code=$? + python3 check.py $ALL_RES "$3" $1 3d 3 || check_code=$? if [[ $check_code != 0 ]]; then echo MHD modes test \"$3\" FAIL: $check_code exit_code=1 @@ -52,7 +52,7 @@ conv_2d() { mv mhdmodes.out0.final.phdf mhd_2d_${1}_${res}_end.phdf done check_code=0 - python check.py $ALL_RES "$3" $1 2d || check_code=$? + python3 check.py $ALL_RES "$3" $1 2d || check_code=$? if [[ $check_code != 0 ]]; then echo MHD modes test \"$3\" FAIL: $check_code exit_code=1 @@ -107,13 +107,13 @@ conv_2d alfven_imex_semi "mhdmodes/nmode=2 driver/type=imex GRMHD/implicit=true conv_2d fast_imex_semi "mhdmodes/nmode=3 driver/type=imex GRMHD/implicit=true b_field/implicit=false" "fast mode 3D, ImEx semi-implicit" # KHARMA driver -conv_2d slow_kharma_ct "mhdmodes/nmode=1 driver/type=kharma b_field/solver=face_ct" "slow mode in 2D, KHARMA driver w/face CT" -conv_2d alfven_kharma_ct "mhdmodes/nmode=2 driver/type=kharma b_field/solver=face_ct" "Alfven mode in 2D, KHARMA driver w/face CT" -conv_2d fast_kharma_ct "mhdmodes/nmode=3 driver/type=kharma b_field/solver=face_ct" "fast mode in 2D, KHARMA driver w/face CT" +conv_2d slow_kharma_ct "mhdmodes/nmode=1 driver/type=kharma b_field/solver=face_ct b_field/ct_scheme=bs99" "slow mode in 2D, KHARMA driver w/face CT" +conv_2d alfven_kharma_ct "mhdmodes/nmode=2 driver/type=kharma b_field/solver=face_ct b_field/ct_scheme=bs99" "Alfven mode in 2D, KHARMA driver w/face CT" +conv_2d fast_kharma_ct "mhdmodes/nmode=3 driver/type=kharma b_field/solver=face_ct b_field/ct_scheme=bs99" "fast mode in 2D, KHARMA driver w/face CT" # ImEx driver -conv_2d slow_imex_ct "mhdmodes/nmode=1 driver/type=imex b_field/solver=face_ct" "slow mode in 2D, ImEx explicit w/face CT" -conv_2d alfven_imex_ct "mhdmodes/nmode=2 driver/type=imex b_field/solver=face_ct" "Alfven mode in 2D, ImEx explicit w/face CT" -conv_2d fast_imex_ct "mhdmodes/nmode=3 driver/type=imex b_field/solver=face_ct" "fast mode in 2D, ImEx explicit w/face CT" +conv_2d slow_imex_ct "mhdmodes/nmode=1 driver/type=imex b_field/solver=face_ct b_field/ct_scheme=bs99" "slow mode in 2D, ImEx explicit w/face CT" +conv_2d alfven_imex_ct "mhdmodes/nmode=2 driver/type=imex b_field/solver=face_ct b_field/ct_scheme=bs99" "Alfven mode in 2D, ImEx explicit w/face CT" +conv_2d fast_imex_ct "mhdmodes/nmode=3 driver/type=imex b_field/solver=face_ct b_field/ct_scheme=bs99" "fast mode in 2D, ImEx explicit w/face CT" # Upwinded fluxes conv_2d slow_kharma_ct_gs05_0 "mhdmodes/nmode=1 driver/type=kharma b_field/solver=face_ct b_field/ct_scheme=gs05_0" "slow mode in 2D, KHARMA driver w/epsilon_0 flux" conv_2d alfven_kharma_ct_gs05_0 "mhdmodes/nmode=2 driver/type=kharma b_field/solver=face_ct b_field/ct_scheme=gs05_0" "Alfven mode in 2D, KHARMA driver w/epsilon_0 flux"