Skip to content

Commit

Permalink
Sync primitive variables if they're fundamental
Browse files Browse the repository at this point in the history 
A few commits back I removed sync_prims, reasoning conserved variables
should be prolongated/restricted, and could always be recovered from
each other.  Both points were wrong: primitive vars can be prolongated
and restricted on boundaries fine, though the latter is not ideal. And,
in EMHD, it is not straightforward to recover P from U, as this is
normally done inline with the computation of the next step's state.

This commit switches to syncing primitive variables (sync_prims) anytime
they're fundamental (ImEx and simple drivers, "prims_are_fundamental")
Note the "conserved" B field is *always* what is sync'd, regardless
of the other primitive or conserved variables.

It also avoids loading the inverter package if GRMHD is implicitly
evolved, and expands some computed domains related to B to work
at the last prolongation operator bug before AMR.
  • Loading branch information
Ben Prather committed Oct 10, 2023
1 parent 20c2a11 commit 0541185
Show file tree
Hide file tree
Showing 10 changed files with 96 additions and 131 deletions.
36 changes: 20 additions & 16 deletions kharma/b_ct/b_ct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
using namespace parthenon;
using parthenon::refinement_ops::ProlongateSharedMinMod;
using parthenon::refinement_ops::RestrictAverage;
using parthenon::refinement_ops::ProlongateInternalAverage;

std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages)
{
Expand Down Expand Up @@ -93,6 +94,8 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
m = Metadata(flags_cons_f);
if (!lazy_prolongation)
m.RegisterRefinementOps<ProlongateSharedMinMod, RestrictAverage, ProlongateInternalOlivares>();
else
m.RegisterRefinementOps<ProlongateSharedMinMod, RestrictAverage, ProlongateInternalAverage>();
pkg->AddField("cons.fB", m);

// Cell-centered versions. Needed for BS, not for other schemes.
Expand Down Expand Up @@ -164,7 +167,7 @@ TaskStatus B_CT::MeshUtoP(MeshData<Real> *md, IndexDomain domain, bool coarse)
return TaskStatus::complete;
}

void B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
TaskStatus B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
{
auto pmb = rc->GetBlockPointer();
const int ndim = pmb->pmy_mesh->ndim;
Expand Down Expand Up @@ -204,6 +207,8 @@ void B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
B_U(v, k, j, i) = B_P(v, k, j, i) * G.gdet(Loci::center, j, i);
}
);

return TaskStatus::complete;
}

TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
Expand All @@ -215,12 +220,9 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
auto& emf_pack = md->PackVariables(std::vector<std::string>{"B_CT.emf"});

// Figure out indices
const IndexRange3 b = KDomain::GetRange(md, IndexDomain::interior, 0, 0);
const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, -1, 2);
const IndexRange3 b = KDomain::GetRange(md, IndexDomain::entire, 0, 0);
const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::entire, 1, 1);
const IndexRange block = IndexRange{0, emf_pack.GetDim(5)-1};
const int kd = ndim > 2 ? 1 : 0;
const int jd = ndim > 1 ? 1 : 0;
const int id = ndim > 0 ? 1 : 0;

auto pmb0 = md->GetBlockData(0)->GetBlockPointer().get();

Expand Down Expand Up @@ -269,28 +271,30 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
auto& B_U = md->PackVariablesAndFluxes(std::vector<std::string>{"cons.B"});
auto& B_P = md->PackVariables(std::vector<std::string>{"prims.B"});
// emf in center == -v x B
const IndexRange3 bc = KDomain::GetRange(md, IndexDomain::entire);
pmb0->par_for("B_CT_emfc", block.s, block.e, bc.ks, bc.ke, bc.js, bc.je, bc.is, bc.ie,
pmb0->par_for("B_CT_emfc", block.s, block.e, b.ks, b.ke, b.js, b.je, b.is, b.ie,
KOKKOS_LAMBDA (const int &bl, const int &k, const int &j, const int &i) {
VLOOP emfc(bl, v, k, j, i) = 0.;
VLOOP3 emfc(bl, x, k, j, i) -= antisym(v, w, x) * uvec(bl, v, k, j, i) * B_U(bl, w, k, j, i);
}
);

if (scheme == "gs05_0") {
const int kd = ndim > 2 ? 1 : 0;
const int jd = ndim > 1 ? 1 : 0;
const int id = ndim > 0 ? 1 : 0;
pmb0->par_for("B_CT_emf_GS05_0", block.s, block.e, b1.ks, b1.ke, b1.js, b1.je, b1.is, b1.ie,
KOKKOS_LAMBDA (const int &bl, const int &k, const int &j, const int &i) {
const auto& G = B_U.GetCoords(bl);
// Just subtract centered emf from twice the face version
// More stable for planar flows even without anything fancy
emf_pack(bl, E1, 0, k, j, i) = 2 * emf_pack(bl, E1, 0, k, j, i)
- 0.25*(emfc(bl, V1, k, j, i) + emfc(bl, V1, k, j - jd, i)
+ emfc(bl, V1, k, j - jd, i) + emfc(bl, V1, k - kd, j - jd, i));
- 0.25*(emfc(bl, V1, k, j, i) + emfc(bl, V1, k, j - jd, i)
+ emfc(bl, V1, k, j - jd, i) + emfc(bl, V1, k - kd, j - jd, i));
emf_pack(bl, E2, 0, k, j, i) = 2 * emf_pack(bl, E2, 0, k, j, i)
- 0.25*(emfc(bl, V2, k, j, i) + emfc(bl, V2, k, j, i - id)
+ emfc(bl, V2, k - kd, j, i) + emfc(bl, V2, k - kd, j, i - id));
- 0.25*(emfc(bl, V2, k, j, i) + emfc(bl, V2, k, j, i - id)
+ emfc(bl, V2, k - kd, j, i) + emfc(bl, V2, k - kd, j, i - id));
emf_pack(bl, E3, 0, k, j, i) = 2 * emf_pack(bl, E3, 0, k, j, i)
- 0.25*(emfc(bl, V3, k, j, i) + emfc(bl, V3, k, j, i - id)
- 0.25*(emfc(bl, V3, k, j, i) + emfc(bl, V3, k, j, i - id)
+ emfc(bl, V3, k, j - jd, i) + emfc(bl, V3, k, j - jd, i - id));
}
);
Expand All @@ -301,14 +305,14 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
pmb0->par_for("B_CT_emf_GS05_c", block.s, block.e, b1.ks, b1.ke, b1.js, b1.je, b1.is, b1.ie,
KOKKOS_LAMBDA (const int &bl, const int &k, const int &j, const int &i) {
const auto& G = B_U.GetCoords(bl);

// "simple" flux + upwinding method, Stone & Gardiner '09 but also in Stone+08 etc.
// Upwinded differences take in order (1-indexed):
// 1. EMF component direction to calculate
// 2. Direction of derivative
// 3. Direction of upwinding
// ...then zone number...
// and finally, a boolean indicating a leftward (e.g., i-3/4) vs rightward (i-1/4) position
// TODO(BSP) This doesn't properly support 2D. Yell when it's chosen?
if (ndim > 2) {
emf_pack(bl, E1, 0, k, j, i) +=
0.25*(upwind_diff(B_U(bl), emfc(bl), uvecf(bl), 1, 3, 2, k, j, i, false)
Expand Down Expand Up @@ -344,8 +348,8 @@ TaskStatus B_CT::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt)
auto& emf_pack = md->PackVariables(std::vector<std::string>{"B_CT.emf"});

// Figure out indices
const IndexRange3 b = KDomain::GetRange(md, IndexDomain::interior, 0, 0);
const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, 0, 1);
const IndexRange3 b = KDomain::GetRange(md, IndexDomain::entire, 0, 0);
const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::entire, 0, 1);
const IndexRange block = IndexRange{0, emf_pack.GetDim(5)-1};

auto pmb0 = md->GetBlockData(0)->GetBlockPointer().get();
Expand Down
14 changes: 5 additions & 9 deletions kharma/b_ct/b_ct.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
* input: Conserved B = sqrt(-gdet) * B^i
* output: Primitive B = B^i
*/
void BlockUtoP(MeshBlockData<Real> *mbd, IndexDomain domain, bool coarse=false);
TaskStatus BlockUtoP(MeshBlockData<Real> *mbd, IndexDomain domain, bool coarse=false);
TaskStatus MeshUtoP(MeshData<Real> *md, IndexDomain domain, bool coarse=false);

/**
Expand Down Expand Up @@ -283,9 +283,10 @@ struct ProlongateInternalOlivares {
const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s;

// Coefficients selecting a particular formula (see Olivares et al. 2019)
// TODO options here. This corresponds to Cunningham, but we could have:
// 1. differences of squares of zone dimesnions (Toth)
// 2. heuristic based on flux difference of top vs bottom halves (Olivares)
// TODO options here. There are 3 presented:
// 1. Zeros (Cunningham)
// 2. differences of squares of zone dimesnions (Toth)
// 3. heuristic based on flux difference of top vs bottom halves (Olivares)
// constexpr Real a[3] = {0., 0., 0.};
const Real a[3] = {(SQR(coords.Dxc<2>(fj)) - SQR(coords.Dxc<3>(fk))) / (SQR(coords.Dxc<2>(fj)) + SQR(coords.Dxc<3>(fk))),
(SQR(coords.Dxc<3>(fk)) - SQR(coords.Dxc<1>(fi))) / (SQR(coords.Dxc<3>(fk)) + SQR(coords.Dxc<1>(fi))),
Expand Down Expand Up @@ -322,11 +323,6 @@ struct ProlongateInternalOlivares {
+ coeff[elem][2]*F<third,me,-1,DIM>(fine, coords, l, m, n, fk, fj, fi)
+ coeff[elem][3]*F<third,me,next,DIM>(fine, coords, l, m, n, fk, fj, fi))
) / coords.Volume<el>(fk+off_k, fj+off_j, fi+off_i);
//printf("%d %d\n", fi, fj);
// if (fi == 56 && fj == 70)
// printf("I used dir %d offset %d %d %d, %d %d %d\n", me+1,
// off_k-diff_k, off_j-diff_j, off_i-diff_i,
// off_k+diff_k, off_j+diff_j, off_i+diff_i);
}
}
};
Expand Down
24 changes: 7 additions & 17 deletions kharma/b_flux_ct/b_flux_ct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,22 +101,11 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
: Metadata::GetUserFlag("Explicit");

// Flags for B fields
std::vector<MetadataFlag> flags_b = {Metadata::Cell, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};

// "primitive" B field is field, "conserved" is flux
auto flags_prim = packages->Get("Driver")->Param<std::vector<MetadataFlag>>("prim_flags");
flags_prim.insert(flags_prim.end(), flags_b.begin(), flags_b.end());
auto flags_cons = packages->Get("Driver")->Param<std::vector<MetadataFlag>>("cons_flags");
flags_cons.insert(flags_cons.end(), flags_b.begin(), flags_b.end());

// Always sync B field conserved var, for standardization with B_CT
// god std::vector is verbose
if (std::find(flags_cons.begin(), flags_cons.end(), Metadata::FillGhost) == flags_cons.end()) {
flags_cons.push_back(Metadata::FillGhost);
}
if (std::find(flags_prim.begin(), flags_prim.end(), Metadata::FillGhost) != flags_prim.end()) {
flags_prim.erase(std::remove(flags_prim.begin(), flags_prim.end(), Metadata::FillGhost), flags_prim.end());
}
// We always mark conserved B to be sync'd for consistency, since it's strictly required for B_CT/AMR
std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Derived, Metadata::GetUserFlag("Primitive"),
Metadata::Cell, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Independent, Metadata::Restart, Metadata::FillGhost, Metadata::WithFluxes, Metadata::Conserved,
Metadata::Cell, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};

auto m = Metadata(flags_prim, s_vector);
pkg->AddField("prims.B", m);
Expand Down Expand Up @@ -192,7 +181,7 @@ void MeshUtoP(MeshData<Real> *md, IndexDomain domain, bool coarse)
}
);
}
void BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
TaskStatus BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
{
auto pmb = rc->GetBlockPointer();

Expand All @@ -213,6 +202,7 @@ void BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
B_P(mu, k, j, i) = B_U(mu, k, j, i) / G.gdet(Loci::center, j, i);
}
);
return TaskStatus::complete;
}

void MeshPtoU(MeshData<Real> *md, IndexDomain domain, bool coarse)
Expand Down
2 changes: 1 addition & 1 deletion kharma/b_flux_ct/b_flux_ct.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
* input: Conserved B = sqrt(-gdet) * B^i
* output: Primitive B = B^i
*/
void BlockUtoP(MeshBlockData<Real> *md, IndexDomain domain, bool coarse=false);
TaskStatus BlockUtoP(MeshBlockData<Real> *md, IndexDomain domain, bool coarse=false);
void MeshUtoP(MeshData<Real> *md, IndexDomain domain, bool coarse=false);

/**
Expand Down
52 changes: 27 additions & 25 deletions kharma/boundaries/boundaries.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,12 +36,15 @@
#include "decs.hpp"
#include "domain.hpp"
#include "kharma.hpp"
#include "flux.hpp"
#include "flux_functions.hpp"
#include "grmhd_functions.hpp"
#include "pack.hpp"
#include "types.hpp"

#include "b_ct.hpp"
#include "b_flux_ct.hpp"
#include "flux.hpp"

// Parthenon's boundaries
#include <bvals/boundary_conditions.hpp>

Expand Down Expand Up @@ -258,13 +261,6 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
const auto btype_name = params.Get<std::string>(bname);
const auto bdir = BoundaryDirection(bface);

// If we're pretending to sync primitives, but applying physical bounds
// to conserved variables, make sure we're up to date
if (pmb->packages.Get<KHARMAPackage>("Driver")->Param<bool>("prims_are_fundamental") &&
params.Get<bool>("domain_bounds_on_conserved")) {
Flux::BlockPtoU_Send(rc.get(), domain, coarse);
}

Flag("Apply "+bname+" boundary: "+btype_name);
pkg->KBoundaries[bface](rc, coarse);
EndFlag();
Expand Down Expand Up @@ -340,25 +336,31 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
}
}

// If we applied the domain boundary to primitives (as we usually do)...
if (!params.Get<bool>("domain_bounds_on_conserved")) {
bool sync_prims = rc->GetBlockPointer()->packages.Get("Driver")->Param<bool>("sync_prims");
// There are two modes of operation here:
if (sync_prims) {
// 1. ImEx w/o AMR:
// PRIMITIVE variables (only) are marked FillGhost
// So, run PtoU on EVERYTHING (and correct the B field)
CorrectBPrimitive(rc, domain, coarse);
Flux::BlockPtoU(rc.get(), domain, coarse);
} else {
// 2. Normal (KHARMA driver, ImEx w/AMR):
// CONSERVED variables are marked FillGhost, plus FLUID PRIMITIVES.
// So, run PtoU on FLUID, and UtoP on EVERYTHING ELSE
Packages::BoundaryPtoUElseUtoP(rc.get(), domain, coarse);
bool sync_prims = rc->GetBlockPointer()->packages.Get("Driver")->Param<bool>("sync_prims");
// There are two modes of operation here:
if (sync_prims) {
// 1. Exchange/prolongate/restrict PRIMITIVE variables: (ImEx driver)
// Primitive variables and conserved B field are marked FillGhost
// Explicitly run UtoP on B field, then PtoU on everything
// TODO there should be a set of B field wrappers that dispatch this
auto pkgs = pmb->packages.AllPackages();
if (pkgs.count("B_FluxCT")) {
B_FluxCT::BlockUtoP(rc.get(), IndexDomain::entire);
} else if (pkgs.count("B_CT")) {
B_CT::BlockUtoP(rc.get(), IndexDomain::entire);
}
Flux::BlockPtoU(rc.get(), domain, coarse);
} else {
// These get applied the same way regardless of driver
Packages::BoundaryUtoP(rc.get(), domain, coarse);
// 2. Exchange/prolongate/restrict CONSERVED variables: (KHARMA driver, maybe ImEx+AMR)
// Conserved variables are marked FillGhost, plus FLUID PRIMITIVES.
if (!params.Get<bool>("domain_bounds_on_conserved")) {
// To apply primitive boundaries to GRMHD, we run PtoU on that ONLY,
// and UtoP on EVERYTHING ELSE
Packages::BoundaryPtoUElseUtoP(rc.get(), domain, coarse);
} else {
// If we want to apply boundaries to conserved vars, just run UtoP on EVERYTHING
Packages::BoundaryUtoP(rc.get(), domain, coarse);
}
}

EndFlag();
Expand Down
Loading

0 comments on commit 0541185

Please sign in to comment.