Skip to content

Commit

Permalink
Merge pull request #91 from AFD-Illinois/dev
Browse files Browse the repository at this point in the history 
KHARMA 2024.5
  • Loading branch information
@bprather
bprather authored May 30, 2024
2 parents 26b8bc8 + 8c5b886 commit 76ae843
Show file tree
Hide file tree
Showing 63 changed files with 1,384 additions and 666 deletions.
2 changes: 0 additions & 2 deletions external/kokkos-kernels/KokkosBatched_Dot_Internal.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,6 @@

#include "KokkosBatched_Util.hpp"

#define KOKKOS_IMPL_DO_NOT_USE_PRINTF(...) ::printf(__VA_ARGS__)

namespace KokkosBatched {

///
Expand Down
13 changes: 10 additions & 3 deletions kharma/b_cd/b_cd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages)
{
throw std::runtime_error("Constraint-damping transport is not functional with modern B field initialization!");

auto pkg = std::make_shared<KHARMAPackage>("B_CD");
Params &params = pkg->AllParams();

Expand Down Expand Up @@ -99,7 +101,12 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
// add callbacks for HST output to the Params struct, identified by the `hist_param_key`
pkg->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<Real> *rc, IndexDomain domain, bool coarse)
Expand Down
4 changes: 4 additions & 0 deletions kharma/b_cd/b_cd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 5 additions & 3 deletions kharma/b_cleanup/b_cleanup.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,11 +139,13 @@ std::shared_ptr<KHARMAPackage> 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);
Expand Down
135 changes: 3 additions & 132 deletions kharma/b_ct/b_ct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,6 @@
#include "grmhd.hpp"
#include "grmhd_functions.hpp"
#include "kharma.hpp"
// TODO eliminate sync
#include "kharma_driver.hpp"

#include <parthenon/parthenon.hpp>
#include <prolong_restrict/pr_ops.hpp>
Expand Down Expand Up @@ -71,7 +69,7 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared

// TODO gs05_alpha, LDZ04 UCT1, LDZ07 UCT2
std::vector<std::string> 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
Expand Down Expand Up @@ -445,132 +443,6 @@ TaskStatus B_CT::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt, IndexDomai
return TaskStatus::complete;
}

void B_CT::ZeroEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &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<Real> *rc, IndexDomain domain, const VariablePack<Real> &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<double> 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<Real> *rc, IndexDomain domain, const VariablePack<Real> &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<F1>(k, j, last_rank_f)
- (fpack(F2, 0, k, j + 1, last_rank_c) * G.Volume<F2>(k, j + 1, last_rank_c)
- fpack(F2, 0, k, j, last_rank_c) * G.Volume<F2>(k, j, last_rank_c));
if (ndim > 2)
new_face -= fpack(F3, 0, k + 1, j, last_rank_c) * G.Volume<F3>(k + 1, j, last_rank_c)
- fpack(F3, 0, k, j, last_rank_c) * G.Volume<F3>(k, j, last_rank_c);

fpack(F1, 0, k, j, i) = outward_sign * new_face / G.Volume<F1>(k, j, i);
}
);
}
} else {
throw std::runtime_error("Divergence-free outflow replacement only implemented in X1!");
}
}

double B_CT::MaxDivB(MeshData<Real> *md)
{
auto pmesh = md->GetMeshPointer();
Expand All @@ -591,7 +463,7 @@ double B_CT::MaxDivB(MeshData<Real> *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);
Expand All @@ -614,14 +486,13 @@ double B_CT::BlockMaxDivB(MeshBlockData<Real> *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<Real> *md, bool all_reduce)
{
if (all_reduce) {
Expand Down
28 changes: 26 additions & 2 deletions kharma/b_ct/b_ct.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,17 +122,41 @@ void CalcDivB(MeshData<Real> *md, std::string divb_field_name="divB");
*
* *mostly. I think
*/
void ZeroEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &emfpack, bool coarse);
void ZeroBoundaryEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &emfpack, bool coarse);

/**
* Average all EMFs corresponding to the coordinate pole location, e.g. usually all E1 on X2 faces
*/
void AverageEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &emfpack, bool coarse);
void AverageBoundaryEMF(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &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<Real> *rc, IndexDomain domain, const VariablePack<Real> &emfpack, bool coarse);

/**
* Reset an outflow condition to have no divergence, even if a field line exits the domain.
* Could maybe be used on other boundaries, but resets the perpendicular face so use with caution.
*/
void DestructiveBoundaryClean(MeshBlockData<Real> *rc, IndexDomain domain, const VariablePack<Real> &fpack, bool coarse);

/**
* Take the curl over the whole domain. Defined in-header since it's templated on face and NDIM
*/
template<TE el, int NDIM>
inline void EdgeCurl(MeshBlockData<Real> *rc, const GridVector& A,
const VariablePack<Real>& 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<el, NDIM>(G, A, B_U, k, j, i);
}
);
}

}
Loading

0 comments on commit 76ae843

Please sign in to comment.