diff --git a/kharma/b_ct/b_ct.cpp b/kharma/b_ct/b_ct.cpp index 95fb1414..b4469afb 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); @@ -168,7 +175,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/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/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..79241db3 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); } } @@ -112,7 +114,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 +127,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")) { 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; }