From 6f011073fb53fa7db2c8ca3313f206c70dc016ef Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Fri, 7 Feb 2025 13:37:14 -0700 Subject: [PATCH 1/2] Add an option to read primitive cell-centered B when restarting, for restarting from .phdf files. --- kharma/b_ct/b_ct.cpp | 7 +++++++ kharma/b_flux_ct/b_flux_ct.cpp | 6 ++++++ kharma/prob/post_initialize.cpp | 6 +++++- 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/kharma/b_ct/b_ct.cpp b/kharma/b_ct/b_ct.cpp index 95fb1414..e9a92e82 100644 --- a/kharma/b_ct/b_ct.cpp +++ b/kharma/b_ct/b_ct.cpp @@ -94,6 +94,13 @@ std::shared_ptr B_CT::Initialize(ParameterInput *pin, std::shared Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector}; std::vector flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::Conserved, Metadata::WithFluxes, Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector}; + + // KHARMA now (vaguely) supports restarting from dump (.phdf) files. + // To do so we need to read cell-centered primitive B and interpolate to faces, as we do with iharm3d restart files + if (pin->GetOrAddBoolean("b_field", "restart_from_prims", false)) { + flags_prim.push_back(Metadata::Restart); + } + std::vector s_vector({NVEC}); m = Metadata(flags_prim, s_vector); pkg->AddField("prims.B", m); diff --git a/kharma/b_flux_ct/b_flux_ct.cpp b/kharma/b_flux_ct/b_flux_ct.cpp index 9f410a54..91bd2e98 100644 --- a/kharma/b_flux_ct/b_flux_ct.cpp +++ b/kharma/b_flux_ct/b_flux_ct.cpp @@ -112,6 +112,12 @@ std::shared_ptr Initialize(ParameterInput *pin, std::shared_ptr

flags_cons = {Metadata::Real, Metadata::Independent, Metadata::Restart, Metadata::FillGhost, Metadata::WithFluxes, Metadata::Conserved, Metadata::Cell, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector}; + // KHARMA now (vaguely) supports restarting from dump (.phdf) files. + // To do so we need to read cell-centered primitive B and interpolate to faces, as we do with iharm3d restart files + if (pin->GetOrAddBoolean("b_field", "restart_from_prims", false)) { + flags_prim.push_back(Metadata::Restart); + } + auto m = Metadata(flags_prim, s_vector); pkg->AddField("prims.B", m); m = Metadata(flags_cons, s_vector); diff --git a/kharma/prob/post_initialize.cpp b/kharma/prob/post_initialize.cpp index 2c381e63..65f3235e 100644 --- a/kharma/prob/post_initialize.cpp +++ b/kharma/prob/post_initialize.cpp @@ -112,7 +112,10 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) // We only record the conserved magnetic field in KHARMA restarts, // but we record primitive field in iharm3d restarts bool iharm3d_restart = prob_name == "resize_restart"; - if (!iharm3d_restart) { + // but sometimes in a pinch we want to restart from .phdf files, which are like iharm3d restarts + bool prims_only_restart = pin->GetBoolean("b_field", "restart_from_prims"); + if (!(iharm3d_restart || prims_only_restart)) { + std::cerr << "Restoring B field from conserved values" << std::endl; if (pkgs.count("B_FluxCT")) { B_FluxCT::MeshUtoP(md.get(), IndexDomain::entire); } else if (pkgs.count("B_CT")) { @@ -122,6 +125,7 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) EMHD::MeshUtoP(md.get(), IndexDomain::entire); } } else { + std::cerr << "Restoring B field from primitive values" << std::endl; if (pkgs.count("B_FluxCT")) { B_FluxCT::MeshPtoU(md.get(), IndexDomain::entire); } else if (pkgs.count("B_CT")) { From 6e2dc9c32ed54a6221c13f0e6dd549faf067bb1d Mon Sep 17 00:00:00 2001 From: Ben Prather Date: Wed, 19 Feb 2025 16:09:59 -0700 Subject: [PATCH 2/2] Update first interior cells on B_CT::UtoP Previously, B_CT::UtoP on e.g. IndexDomain::inner_x1 would only update ghost cell centers. But changing the flux on the domain face calls for updating the first interior row too, since it's involved in averaging. Also patches up initializing/normalizing B in the presence of boundaries --- kharma/b_ct/b_ct.cpp | 12 +++++++++++- kharma/boundaries/boundaries.cpp | 2 ++ kharma/driver/kharma_driver.cpp | 24 ++++++++++++------------ kharma/prob/post_initialize.cpp | 2 ++ kharma/prob/seed_B.cpp | 31 +++++++++++++++---------------- 5 files changed, 42 insertions(+), 29 deletions(-) diff --git a/kharma/b_ct/b_ct.cpp b/kharma/b_ct/b_ct.cpp index 95fb1414..297b83ee 100644 --- a/kharma/b_ct/b_ct.cpp +++ b/kharma/b_ct/b_ct.cpp @@ -168,7 +168,17 @@ TaskStatus B_CT::BlockUtoP(MeshBlockData *rc, IndexDomain domain, bool coa // Return if we're not syncing U & P at all (e.g. edges) if (B_Uf.GetDim(4) == 0) return TaskStatus::complete; - const IndexRange3 bc = KDomain::GetRange(rc, domain, coarse); + // We need one cell inside of the domain, since it will be updated by the domain-face field + // this oversteps "interior" by a zone + IndexRange3 bct; + if (domain == IndexDomain::interior || domain == IndexDomain::entire) { + bct = KDomain::GetRange(rc, domain, coarse); + } else { + const BoundaryFace bface = KBoundaries::BoundaryFaceOf(domain); + const bool binner = KBoundaries::BoundaryIsInner(bface); + bct = KDomain::GetRange(rc, domain, (binner) ? 0 : -1, (binner) ? 1 : 0, coarse); + } + const IndexRange3 bc = bct; // Average the primitive vals to zone centers pmb->par_for("UtoP_B_center", bc.ks, bc.ke, bc.js, bc.je, bc.is, bc.ie, diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp index ac289b79..002089eb 100644 --- a/kharma/boundaries/boundaries.cpp +++ b/kharma/boundaries/boundaries.cpp @@ -716,6 +716,8 @@ TaskStatus KBoundaries::FixFlux(MeshData *md) const Real gam = pmb->packages.Get("GRMHD")->Param("gamma"); const auto& G = pmb->coords; const int nvar = F.GetDim(4); + // Loci "outer" and "inner" refer to right and left. + // At right edge, our "outer" toward domain is Loci "inner" as in left half const Loci loc = (binner) ? Loci::outer_half : Loci::inner_half; const IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, CC); diff --git a/kharma/driver/kharma_driver.cpp b/kharma/driver/kharma_driver.cpp index 921098b9..05a6bc5f 100644 --- a/kharma/driver/kharma_driver.cpp +++ b/kharma/driver/kharma_driver.cpp @@ -1,26 +1,26 @@ -/* +/* * File: kharma_driver.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 @@ -85,8 +85,8 @@ std::shared_ptr KHARMADriver::Initialize(ParameterInput *pin, std // 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 + // 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"); @@ -187,14 +187,14 @@ TaskStatus KHARMADriver::SyncAllBounds(std::shared_ptr> &md) Flag("SyncAllBounds"); TaskID t_none(0); - //MPIBarrier(); + MPIBarrier(); TaskCollection tc; auto tr = tc.AddRegion(1); AddBoundarySync(t_none, tr[0], md); while (!tr.Execute()); - //MPIBarrier(); + MPIBarrier(); EndFlag(); return TaskStatus::complete; diff --git a/kharma/prob/post_initialize.cpp b/kharma/prob/post_initialize.cpp index 2c381e63..3caf7226 100644 --- a/kharma/prob/post_initialize.cpp +++ b/kharma/prob/post_initialize.cpp @@ -88,6 +88,8 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart) // normalize the magnetic field according to the max density bool is_torus = prob_name == "torus"; if (pin->GetOrAddBoolean("b_field", "norm", is_torus)) { + // Sync to apply physical boundary conditions and fill B ghosts correctly + KHARMADriver::SyncAllBounds(md); NormalizeBField(md.get(), pin); } } diff --git a/kharma/prob/seed_B.cpp b/kharma/prob/seed_B.cpp index 4b6b4b58..d0b15b86 100644 --- a/kharma/prob/seed_B.cpp +++ b/kharma/prob/seed_B.cpp @@ -1,25 +1,25 @@ -/* +/* * File: seed_B.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 @@ -47,19 +47,20 @@ using namespace parthenon; /** * Perform a Parthenon MPI reduction. - * Should only be used in initialization code, as the - * reducer object & MPI comm are created on entry & - * cleaned on exit + * Should only be used in initialization code, as it + * is very not fast. */ template inline T MPIReduce_once(T f, MPI_Op O) { + MPIBarrier(); parthenon::AllReduce reduction; reduction.val = f; reduction.StartReduce(O); // Wait on results while (reduction.CheckReduce() == parthenon::TaskStatus::incomplete); // TODO catch errors? + MPIBarrier(); return reduction.val; } @@ -114,7 +115,7 @@ TaskStatus SeedBFieldType(MeshBlockData *rc, ParameterInput *pin, IndexDom if constexpr (Seed == BSeedType::constant || Seed == BSeedType::monopole || Seed == BSeedType::orszag_tang || - Seed == BSeedType::wave || + Seed == BSeedType::wave || Seed == BSeedType::shock_tube) { // All custom B fields should set what they need of these. @@ -331,7 +332,8 @@ TaskStatus SeedBFieldType(MeshBlockData *rc, ParameterInput *pin, IndexDom // So, we preserve exact values in the no-tilt case. A(V3, k, j, i) = Aphi; } - }); + } + ); if (pkgs.count("B_CT")) { auto B_Uf = rc->PackVariables(std::vector{"cons.fB"}); @@ -526,9 +528,6 @@ TaskStatus NormalizeBField(MeshData *md, ParameterInput *pin) } } - // We've been initializing/manipulating P - Flux::MeshPtoU(md, IndexDomain::entire); - EndFlag(); //NormBField return TaskStatus::complete; }