Skip to content

Commit

Permalink
Merge branch 'feature/meshify' into feature/parthenon-bump
Browse files Browse the repository at this point in the history
  • Loading branch information
Ben Prather committed Nov 17, 2023
2 parents 1f65b7e + abb6d90 commit 0223f8a
Show file tree
Hide file tree
Showing 22 changed files with 200 additions and 150 deletions.
33 changes: 1 addition & 32 deletions kharma/b_ct/b_ct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,9 +174,6 @@ TaskStatus B_CT::BlockUtoP(MeshBlockData<Real> *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;

// TODO get rid of prims on faces probably

// Update the primitive B-fields on faces
const IndexRange3 bc = KDomain::GetRange(rc, domain, coarse);

// Average the primitive vals to zone centers
Expand Down Expand Up @@ -327,6 +324,7 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
} else {
throw std::invalid_argument("Invalid CT scheme specified! Must be one of bs99, gs05_0, gs05_c!");
}

return TaskStatus::complete;
}

Expand Down Expand Up @@ -380,35 +378,6 @@ TaskStatus B_CT::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt)
}
);

// Explicitly zero polar faces
// In spherical, zero B2 on X2 face regardless of boundary condition
// This shouldn't interfere with divB since the face size is zero anyway
if (mdudt->GetBlockData(0)->GetBlockPointer()->coords.coords.is_spherical()) {
const IndexRange ib = mdudt->GetBoundsI(IndexDomain::entire);
const IndexRange kb = mdudt->GetBoundsK(IndexDomain::entire);
const int js = mdudt->GetBoundsJ(IndexDomain::interior).s;
const int je = mdudt->GetBoundsJ(IndexDomain::interior).e + 1; // Face
for (int i_block = 0; i_block < mdudt->NumBlocks(); i_block++) {
auto &rc = mdudt->GetBlockData(i_block);
auto pmb = rc->GetBlockPointer();
auto& dB_Uf_dt_block = rc->PackVariables(std::vector<std::string>{"cons.fB"});
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::inner_x2)) {
pmb->par_for("B_CT_zero_B2_in", kb.s, kb.e, js, js, ib.s, ib.e,
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
dB_Uf_dt_block(F2, 0, k, j, i) = 0;
}
);
}
if (KBoundaries::IsPhysicalBoundary(pmb, BoundaryFace::outer_x2)) {
pmb->par_for("B_CT_zero_B2_out", kb.s, kb.e, je, je, ib.s, ib.e,
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
dB_Uf_dt_block(F2, 0, k, j, i) = 0;
}
);
}
}
}

return TaskStatus::complete;
}

Expand Down
10 changes: 5 additions & 5 deletions kharma/b_flux_ct/b_flux_ct.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ inline Real ReducePhi5(MeshData<Real> *md)
* TODO likely better templated, as with all ND stuff
*/
template<typename Global>
KOKKOS_INLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& B_U, const int& b,
KOKKOS_FORCEINLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& B_U, const int& b,
const int& k, const int& j, const int& i, const bool& do_3D)
{
const double norm = (do_3D) ? 0.25 : 0.5;
Expand All @@ -170,7 +170,7 @@ KOKKOS_INLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& B
return norm*term1/G.Dxc<1>(i) + norm*term2/G.Dxc<2>(j) + norm*term3/G.Dxc<3>(k);
}
template<typename Global>
KOKKOS_INLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& P, const VarMap& m_p,
KOKKOS_FORCEINLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& P, const VarMap& m_p,
const int& b, const int& k, const int& j, const int& i,
const bool& do_3D)
{
Expand Down Expand Up @@ -200,7 +200,7 @@ KOKKOS_INLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& P
* Note this is forward-difference, while previous def is backward
*/
template<typename Global>
KOKKOS_INLINE_FUNCTION void center_grad(const GRCoordinates& G, const Global& P, const int& b,
KOKKOS_FORCEINLINE_FUNCTION void center_grad(const GRCoordinates& G, const Global& P, const int& b,
const int& k, const int& j, const int& i, const bool& do_3D,
double& B1, double& B2, double& B3)
{
Expand All @@ -227,7 +227,7 @@ KOKKOS_INLINE_FUNCTION void center_grad(const GRCoordinates& G, const Global& P,
B3 = norm*term3/G.Dxc<3>(k);
}

KOKKOS_INLINE_FUNCTION void averaged_curl_3D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
KOKKOS_FORCEINLINE_FUNCTION void averaged_curl_3D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
const int& k, const int& j, const int& i)
{
// Take a flux-ct step from the corner potentials.
Expand Down Expand Up @@ -270,7 +270,7 @@ KOKKOS_INLINE_FUNCTION void averaged_curl_3D(const GRCoordinates& G, const GridV
B_U(V3, k, j, i) = (A2c1f - A2c1b) / G.Dxc<1>(i) - (A1c2f - A1c2b) / G.Dxc<2>(j);
}

KOKKOS_INLINE_FUNCTION void averaged_curl_2D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
KOKKOS_FORCEINLINE_FUNCTION void averaged_curl_2D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
const int& k, const int& j, const int& i)
{
// A3,2 derivative
Expand Down
81 changes: 42 additions & 39 deletions kharma/boundaries/boundaries.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,8 +129,6 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
// This is two separate checks, but default to enabling/disabling together for X1 and not elsewhere
bool check_inflow = pin->GetOrAddBoolean("boundaries", "check_inflow_" + bname, check_inflow_global && bdir == X1DIR);
params.Add("check_inflow_" + bname, check_inflow);
bool check_inflow_flux = pin->GetOrAddBoolean("boundaries", "check_inflow_flux_" + bname, check_inflow);
params.Add("check_inflow_flux_" + bname, check_inflow_flux);

// Ensure fluxes through the zero-size face at the pole are zero
bool zero_flux = pin->GetOrAddBoolean("boundaries", "zero_flux_" + bname, zero_polar_flux && bdir == X2DIR);
Expand Down Expand Up @@ -260,11 +258,52 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
const auto bname = BoundaryName(bface);
const auto btype_name = params.Get<std::string>(bname);
const auto bdir = BoundaryDirection(bface);
const bool binner = BoundaryIsInner(bface);

Flag("Apply "+bname+" boundary: "+btype_name);
pkg->KBoundaries[bface](rc, coarse);
EndFlag();

// If we're syncing EMFs and in spherical, explicitly zero polar faces
// TODO allow any other face?
auto& emfpack = rc->PackVariables(std::vector<std::string>{"B_CT.emf"});
if (bdir == X2DIR && pmb->coords.coords.is_spherical() && emfpack.GetDim(4) > 0) {
Flag("BoundaryEdge_"+bname);
for (TE el : {TE::E1, TE::E3}) {
int off = (binner) ? 1 : -1;
pmb->par_for_bndry(
"zero_EMF", IndexRange{0,0}, domain, el, coarse,
KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) {
emfpack(el, v, k, j + off, i) = 0;
}
);
}
EndFlag();
}

// Zero/invert X2 faces at polar X2 boundary
auto fpack = rc->PackVariables({Metadata::Face, Metadata::FillGhost});
if (bdir == X2DIR && pmb->coords.coords.is_spherical() && fpack.GetDim(4) > 0) {
Flag("BoundaryFace_"+bname);
// Zero face fluxes
auto b = KDomain::GetRange(rc, domain, coarse);
// "domain" is the boundary here
auto jf = (binner) ? b.je + 1 : b.js;
pmb->par_for(
"zero_polar_" + bname, b.ks, b.ke, jf, jf, b.is, b.ie,
KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
fpack(F2, 0, k, j, i) = 0.;
}
);
pmb->par_for_bndry(
"invert_F2_" + bname, IndexRange{0, fpack.GetDim(4)-1}, domain, F2, coarse,
KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) {
fpack(F2, 0, k, j, i) *= -1;
}
);
EndFlag();
}

// This will now be called in 2 places we might not expect,
// where we still may want to control the physical bounds:
// 1. Syncing only the EMF during runs with CT
Expand Down Expand Up @@ -376,14 +415,6 @@ void KBoundaries::CheckInflow(std::shared_ptr<MeshBlockData<Real>> &rc, IndexDom
auto P = GRMHD::PackMHDPrims(rc.get(), prims_map, coarse);
const VarMap m_p(prims_map, false);

// Inflow check
// Iterate over zones w/p=0
// pmb->par_for_bndry(
// "check_inflow", IndexRange{0, 0}, domain, CC, coarse,
// KOKKOS_LAMBDA(const int &p, const int &k, const int &j, const int &i) {
// KBoundaries::check_inflow(G, P, domain, m_p.U1, k, j, i);
// }
// );
const auto bface = BoundaryFace(domain);
const auto bname = BoundaryName(bface);
const bool binner = BoundaryIsInner(bface);
Expand All @@ -397,34 +428,6 @@ void KBoundaries::CheckInflow(std::shared_ptr<MeshBlockData<Real>> &rc, IndexDom
);
}

void KBoundaries::CorrectBPrimitive(std::shared_ptr<MeshBlockData<Real>>& rc, IndexDomain domain, bool coarse)
{
Flag("CorrectBPrimitive");
std::shared_ptr<MeshBlock> pmb = rc->GetBlockPointer();
auto B_P = rc->PackVariables(std::vector<std::string>{"prims.B"});
// Return if no field to correct
if (B_P.GetDim(4) == 0) return;

const auto& G = pmb->coords;

const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
const int dir = BoundaryDirection(domain);
const auto &range = (dir == 1) ? bounds.GetBoundsI(IndexDomain::interior)
: (dir == 2 ? bounds.GetBoundsJ(IndexDomain::interior)
: bounds.GetBoundsK(IndexDomain::interior));
const int ref = BoundaryIsInner(domain) ? range.s : range.e;

pmb->par_for_bndry(
"Correct_B_P", IndexRange{0,NVEC-1}, domain, CC, coarse,
KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) {
B_P(v, k, j, i) *= G.gdet(Loci::center, (dir == 2) ? ref : j, (dir == 1) ? ref : i)
/ G.gdet(Loci::center, j, i);
}
);

EndFlag();
}

TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
{
auto pmesh = md->GetMeshPointer();
Expand Down Expand Up @@ -477,7 +480,7 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
auto &F = rc->PackVariablesAndFluxes({Metadata::WithFluxes}, cons_map);

// If we should check inflow on this face...
if (params.Get<bool>("check_inflow_flux_" + bname)) {
if (params.Get<bool>("check_inflow_" + bname)) {
const int m_rho = cons_map["cons.rho"].first;
// ...and if this face of the block corresponds to a global boundary...
if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
Expand Down
5 changes: 0 additions & 5 deletions kharma/boundaries/boundaries.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,11 +84,6 @@ TaskStatus FixFlux(MeshData<Real> *rc);
*/
void CheckInflow(std::shared_ptr<MeshBlockData<Real>> &rc, IndexDomain domain, bool coarse);

/**
* Correct for geometry when applying primitive B field boundaries
*/
void CorrectBPrimitive(std::shared_ptr<MeshBlockData<Real>>& rc, IndexDomain domain, bool coarse);

/**
* Check for velocity toward the simulation domain in a zone, and eliminate it.
*/
Expand Down
4 changes: 2 additions & 2 deletions kharma/decs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ using GReal = double;
enum class Loci{face1=0, face2, face3, center, corner};

// Return the face location corresponding to the direction 'dir'
KOKKOS_INLINE_FUNCTION Loci loc_of(const int& dir)
KOKKOS_FORCEINLINE_FUNCTION Loci loc_of(const int& dir)
{
switch (dir) {
case 0:
Expand All @@ -119,7 +119,7 @@ KOKKOS_INLINE_FUNCTION Loci loc_of(const int& dir)
return Loci::corner;
}
}
KOKKOS_INLINE_FUNCTION int dir_of(const Loci loc)
KOKKOS_FORCEINLINE_FUNCTION int dir_of(const Loci loc)
{
switch (loc) {
case Loci::center:
Expand Down
Loading

0 comments on commit 0223f8a

Please sign in to comment.