Skip to content

Commit

Permalink
a few more changes to make it compilable
Browse files Browse the repository at this point in the history
  • Loading branch information
Hyerin Cho committed Oct 8, 2024
1 parent 5914c1f commit 1b52648
Show file tree
Hide file tree
Showing 6 changed files with 102 additions and 106 deletions.
10 changes: 9 additions & 1 deletion kharma/driver/imex_step.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include "b_flux_ct.hpp"
#include "electrons.hpp"
#include "inverter.hpp"
#include "ismr.hpp"
#include "grmhd.hpp"
#include "wind.hpp"
// Other headers
Expand Down Expand Up @@ -221,6 +222,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
auto t_copy_guess = t_explicit;
auto t_copy_emhd_vars = t_explicit;
if (use_ideal_guess) {
// This copies only q, dP for which we don't already have a guess
t_copy_emhd_vars = tl.AddTask(t_sources, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::GetUserFlag("EMHDVar")}),
md_sub_step_init.get(), md_solver.get());
t_copy_guess = t_copy_emhd_vars;
Expand Down Expand Up @@ -250,7 +252,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
// If we're entirely explicit, we just declare these equal
auto t_implicit_c = tl.AddTask(t_implicit_step, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::Cell}),
md_solver.get(), md_sub_step_final.get());
t_implicit = tl.AddTask(t_implicit_step, WeightedSumDataFace, std::vector<MetadataFlag>({Metadata::Face}),
t_implicit = tl.AddTask(t_implicit_step, WeightedSumDataFace<MetadataFlag>, std::vector<MetadataFlag>({Metadata::Face}),
md_solver.get(), md_solver.get(), 1.0, 0.0, md_sub_step_final.get());
}

Expand Down Expand Up @@ -306,6 +308,12 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
auto t_ptou = tl.AddTask(t_heat_electrons, Flux::MeshPtoU, md_sub_step_final.get(), IndexDomain::entire, false);

auto t_step_done = t_ptou;
if (pkgs.count("ISMR")) {
auto t_derefine_b = t_ptou;
if (pkgs.count("B_CT"))
t_derefine_b = tl.AddTask(t_ptou, B_CT::DerefinePoles, md_sub_step_final.get());
//t_step_done = tl.AddTask(t_derefine_b, ISMR::DerefinePoles, md_sub_step_final.get()); // TODO: this needs to be updated
}

// Estimate next time step based on ctop
if (stage == integrator->nstages) {
Expand Down
8 changes: 4 additions & 4 deletions kharma/driver/kharma_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,7 +382,7 @@ TaskID KHARMADriver::AddStateUpdate(TaskID& t_start, TaskList& tl, MeshData<Real
md_update);
auto t_avg_data = t_avg_data_c;
if (update_face) {
t_avg_data = tl.AddTask(t_start, WeightedSumDataFace,
t_avg_data = tl.AddTask(t_start, WeightedSumDataFace<MetadataFlag>,
std::vector<MetadataFlag>(flags_face),
md_sub_step_init, md_full_step_init,
integrator->gam0[stage-1], integrator->gam1[stage-1],
Expand All @@ -396,7 +396,7 @@ TaskID KHARMADriver::AddStateUpdate(TaskID& t_start, TaskList& tl, MeshData<Real
md_update);
auto t_update = t_update_c;
if (update_face) {
t_update = tl.AddTask(t_avg_data, WeightedSumDataFace,
t_update = tl.AddTask(t_avg_data, WeightedSumDataFace<MetadataFlag>,
std::vector<MetadataFlag>(flags_face),
md_update, md_flux_src,
1.0, integrator->beta[stage-1] * integrator->dt,
Expand Down Expand Up @@ -432,7 +432,7 @@ TaskID KHARMADriver::AddStateUpdateIdealGuess(TaskID& t_start, TaskList& tl, Mes
md_update);
auto t_avg_data = t_avg_data_c;
if (update_face) {
t_avg_data = tl.AddTask(t_start, WeightedSumDataFace,
t_avg_data = tl.AddTask(t_start, WeightedSumDataFace<MetadataFlag>,
std::vector<MetadataFlag>(flags_face),
md_sub_step_init, md_full_step_init,
integrator->gam0[stage-1], integrator->gam1[stage-1],
Expand All @@ -446,7 +446,7 @@ TaskID KHARMADriver::AddStateUpdateIdealGuess(TaskID& t_start, TaskList& tl, Mes
md_update);
auto t_update = t_update_c;
if (update_face) {
t_update = tl.AddTask(t_avg_data, WeightedSumDataFace,
t_update = tl.AddTask(t_avg_data, WeightedSumDataFace<MetadataFlag>,
std::vector<MetadataFlag>(flags_face),
md_update, md_flux_src,
1.0, integrator->beta[stage-1] * integrator->dt,
Expand Down
47 changes: 35 additions & 12 deletions kharma/reductions/reductions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,28 +51,51 @@ namespace Reductions {
std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages);

/**
* Perform a reduction using operation 'op' over a spherical shell at the given zone, measured from left side of
* innermost block in radius.
* As this only runs on innermost blocks, this is intended for accretion/event horizon
* measurements in black hole simulations.
*/
template<Var var, typename T>
T EHReduction(MeshData<Real> *md, UserHistoryOperation op, int zone);

/**
* Perform a reduction using operation 'op' over a given domain
* Perform a reduction using operation 'op' over a given domain.
* Note startx/stopx are in *embedding* coordinates e.g. Kerr-Schild, not MKS/FMKS
*
* This should be used for all 2D shell sums not around the EH:
* Just set equal min/max, 2D slices are detected
*/
template<Var var, typename T>
T DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const GReal startx[3], const GReal stopx[3], int channel=-1);
template<Var var, typename T>
T DomainReduction(MeshData<Real> *md, UserHistoryOperation op, int channel=-1) {
const GReal startx[3] = {std::numeric_limits<Real>::min(), std::numeric_limits<Real>::min(), std::numeric_limits<Real>::min()};
const GReal stopx[3] = {std::numeric_limits<Real>::max(), std::numeric_limits<Real>::max(), std::numeric_limits<Real>::max()};
const GReal startx[3] = {std::numeric_limits<GReal>::min(), std::numeric_limits<GReal>::min(), std::numeric_limits<GReal>::min()};
const GReal stopx[3] = {std::numeric_limits<GReal>::max(), std::numeric_limits<GReal>::max(), std::numeric_limits<GReal>::max()};
return DomainReduction<var, T>(md, op, startx, stopx, channel);
}
template<Var var, typename T>
T ShellReduction(MeshData<Real> *md, UserHistoryOperation op, GReal r, int channel=-1) {
const GReal startx[3] = {r, std::numeric_limits<GReal>::min(), std::numeric_limits<GReal>::min()};
const GReal stopx[3] = {r, std::numeric_limits<GReal>::max(), std::numeric_limits<GReal>::max()};
return DomainReduction<var, T>(md, op, startx, stopx, channel);
}

// Parthenon doesn't allow taking options, so we define some common reductions
template<Var var>
Real SumAt0(MeshData<Real> *md)
{
const GReal r_in = md->GetMeshPointer()->packages.Get("Reductions")->Param<GReal>("domain_r_in");
return Reductions::ShellReduction<var, Real>(md, UserHistoryOperation::sum, r_in);
}
template<Var var>
Real SumAtEH(MeshData<Real> *md)
{
const GReal r_eh = md->GetMeshPointer()->block_list[0]->coords.coords.get_horizon();
return Reductions::ShellReduction<var, Real>(md, UserHistoryOperation::sum, r_eh);
}
template<Var var>
Real SumAt5M(MeshData<Real> *md)
{
return Reductions::ShellReduction<var, Real>(md, UserHistoryOperation::sum, 5.);
}
template<Var var>
Real Total(MeshData<Real> *md)
{
return Reductions::DomainReduction<var, Real>(md, UserHistoryOperation::sum);
}

/**
* Start reductions with a value you have on hand
*/
Expand Down
104 changes: 19 additions & 85 deletions kharma/reductions/reductions_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,84 +118,10 @@ T Reductions::CheckOnAll(MeshData<Real> *md, int channel)
return reducer.val;
}

#define REDUCE_FUNCTION_CALL G, P(b), m_p, U(b), m_u, cmax(b), cmin(b), emhd_params, gam, k, j, i

template<Reductions::Var var, typename T>
T Reductions::EHReduction(MeshData<Real> *md, UserHistoryOperation op, int zone)
{
Flag("EHReduction");
auto pmesh = md->GetMeshPointer();

const auto& pars = pmesh->packages.Get("GRMHD")->AllParams();
const Real gam = pars.Get<Real>("gamma");
const auto& emhd_params = EMHD::GetEMHDParameters(pmesh->packages);

PackIndexMap prims_map, cons_map;
const auto& P = md->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
const auto& U = md->PackVariablesAndFluxes(std::vector<MetadataFlag>{Metadata::Conserved}, cons_map);
const VarMap m_u(cons_map, true), m_p(prims_map, false);
const auto& cmax = md->PackVariables(std::vector<std::string>{"Flux.cmax"});
const auto& cmin = md->PackVariables(std::vector<std::string>{"Flux.cmin"});

auto pmb0 = md->GetBlockData(0)->GetBlockPointer();
IndexRange ib = pmb0->cellbounds.GetBoundsI(IndexDomain::interior);
IndexRange jb = pmb0->cellbounds.GetBoundsJ(IndexDomain::interior);
IndexRange kb = pmb0->cellbounds.GetBoundsK(IndexDomain::interior);
IndexRange block = IndexRange{0, U.GetDim(5) - 1};

T result(0);
int nb = pmesh->GetNumMeshBlocksThisRank();
for (int iblock=0; iblock < nb; iblock++) {
const auto &pmb = pmesh->block_list[iblock];
// Inner-edge blocks only for speed
if (pmb->boundary_flag[parthenon::BoundaryFace::inner_x1] == BoundaryFlag::user) {
const auto& G = pmb->coords;
T block_result;
switch(op) {
case UserHistoryOperation::sum: {
Kokkos::Sum<T> sum_reducer(block_result);
pmb->par_reduce("accretion_sum", iblock, iblock, kb.s, kb.e, jb.s, jb.e, ib.s+zone, ib.s+zone,
KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
local_result += reduction_var<var>(REDUCE_FUNCTION_CALL) * G.Dxc<3>(k) * G.Dxc<2>(j);
}
, sum_reducer);
result += block_result;
break;
}
case UserHistoryOperation::max: {
Kokkos::Max<T> max_reducer(block_result);
pmb->par_reduce("accretion_sum", iblock, iblock, kb.s, kb.e, jb.s, jb.e, ib.s+zone, ib.s+zone,
KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
const T val = reduction_var<var>(REDUCE_FUNCTION_CALL) * G.Dxc<3>(k) * G.Dxc<2>(j);
if (val > local_result) local_result = val;
}
, max_reducer);
if (block_result > result) result = block_result;
break;
}
case UserHistoryOperation::min: {
Kokkos::Min<T> min_reducer(block_result);
pmb->par_reduce("accretion_sum", iblock, iblock, kb.s, kb.e, jb.s, jb.e, ib.s+zone, ib.s+zone,
KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i, T &local_result) {
const T val = reduction_var<var>(REDUCE_FUNCTION_CALL) * G.Dxc<3>(k) * G.Dxc<2>(j);
if (val < local_result) local_result = val;
}
, min_reducer);
if (block_result < result) result = block_result;
break;
}
}
}
}

EndFlag();
return result;
}

#define INSIDE (x[1] > startx1 && x[2] > startx2 && x[3] > startx3) && \
(trivial1 ? x[1] < startx1 + G.Dxc<1>(i) : x[1] < stopx1) && \
(trivial2 ? x[2] < startx2 + G.Dxc<2>(j) : x[2] < stopx2) && \
(trivial3 ? x[3] < startx3 + G.Dxc<3>(k) : x[3] < stopx3)
(trivial1 ? xin[1] < startx1 : x[1] < stopx1) && \
(trivial2 ? xin[2] < startx2 : x[2] < stopx2) && \
(trivial3 ? xin[3] < startx3 : x[3] < stopx3)

// TODO additionally template on return type to avoid counting flags with Reals
template<Reductions::Var var, typename T>
Expand All @@ -210,8 +136,8 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const

// Just pass in everything we might want. Probably slow?
PackIndexMap prims_map, cons_map;
const auto& P = md->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
const auto& U = md->PackVariablesAndFluxes(std::vector<MetadataFlag>{Metadata::Conserved}, cons_map);
const auto& P = md->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map);
const auto& U = md->PackVariablesAndFluxes(std::vector<MetadataFlag>{Metadata::Conserved, Metadata::Cell}, cons_map);
const VarMap m_u(cons_map, true), m_p(prims_map, false);
const auto& cmax = md->PackVariables(std::vector<std::string>{"Flux.cmax"});
const auto& cmin = md->PackVariables(std::vector<std::string>{"Flux.cmin"});
Expand Down Expand Up @@ -243,14 +169,17 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
switch(op) {
case UserHistoryOperation::sum: {
Kokkos::Sum<T> sum_reducer(result);
// TODO these should be nested parallel loops to skip blocks fully outside domain
pmb0->par_reduce("domain_sum", 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, T &local_result) {
const auto& G = U.GetCoords(b);
GReal x[4];
GReal x[GR_DIM], xin[GR_DIM];
G.coord_embed(k, j, i, Loci::center, x);
if (trivial1 || trivial2 || trivial3)
G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
if(INSIDE) {
local_result += reduction_var<var>(REDUCE_FUNCTION_CALL) *
(!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
}
}
, sum_reducer);
Expand All @@ -262,11 +191,13 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
pmb0->par_reduce("domain_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, T &local_result) {
const auto& G = U.GetCoords(b);
GReal x[4];
GReal x[GR_DIM], xin[GR_DIM];
G.coord_embed(k, j, i, Loci::center, x);
if (trivial1 || trivial2 || trivial3)
G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
if(INSIDE) {
const Real val = reduction_var<var>(REDUCE_FUNCTION_CALL) *
(!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
if (val > local_result) local_result = val;
}
}
Expand All @@ -279,11 +210,13 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const
pmb0->par_reduce("domain_min", 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, T &local_result) {
const auto& G = U.GetCoords(b);
GReal x[4];
GReal x[GR_DIM], xin[GR_DIM];
G.coord_embed(k, j, i, Loci::center, x);
if (trivial1 || trivial2 || trivial3)
G.coord_embed(k - trivial3, j - trivial2, i - trivial1, Loci::center, xin);
if(INSIDE) {
const Real val = reduction_var<var>(REDUCE_FUNCTION_CALL) *
(!trivial3) * G.Dxc<3>(k) * (!trivial2) * G.Dxc<2>(j) * (!trivial1) * G.Dxc<1>(i);
((trivial3) ? 1. : G.Dxc<3>(k)) * ((trivial2) ? 1. : G.Dxc<2>(j)) * ((trivial1) ? 1. : G.Dxc<1>(i));
if (val < local_result) local_result = val;
}
}
Expand All @@ -304,3 +237,4 @@ T Reductions::DomainReduction(MeshData<Real> *md, UserHistoryOperation op, const

#undef INSIDE
#undef REDUCE_FUNCTION_CALL
#undef REDUCE_FUNCTION_ARGS
37 changes: 34 additions & 3 deletions kharma/reductions/reductions_variables.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,20 +39,23 @@
#include "emhd.hpp"
#include "flux_functions.hpp"

using namespace parthenon;

// Arguments to computing any variable defined below
#define REDUCE_FUNCTION_ARGS const GRCoordinates& G, const VariablePack<Real>& P, const VarMap& m_p, \
const VariableFluxPack<Real>& U, const VarMap& m_u, \
const VariablePack<Real>& cmax, const VariablePack<Real>& cmin,\
const EMHD::EMHD_parameters& emhd_params, const Real& gam, const int& k, const int& j, const int& i
// Call for passing a particular block's values
#define REDUCE_FUNCTION_CALL G, P(b), m_p, U(b), m_u, cmax(b), cmin(b), emhd_params, gam, k, j, i

using namespace parthenon;

namespace Reductions {

// Add any new reduction variables to this list, then implementations below
// Not elegant, but fast & portable.
// HIPCC doesn't like passing function pointers as we used to do,
// and it doesn't vectorize anyway. Look forward to more of this pattern in the code
enum class Var{phi, bsq, gas_pressure, mag_pressure, beta,
enum class Var{phi, bsq, gas_pressure, beta, rhou0, mix_T00, mix_T01, mix_T02, mix_T03,
mdot, edot, ldot, mdot_flux, edot_flux, ldot_flux, eht_lum, jet_lum,
nan_ctop, zero_ctop, neg_rho, neg_u, neg_rhout};

Expand All @@ -68,6 +71,7 @@ KOKKOS_INLINE_FUNCTION Real reduction_var<Var::phi>(REDUCE_FUNCTION_ARGS)
return 0.5 * m::abs(U(m_u.B1, k, j, i)); // factor of gdet already in cons.B
}

// Basic variables
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::bsq>(REDUCE_FUNCTION_ARGS)
{
Expand All @@ -88,6 +92,33 @@ KOKKOS_INLINE_FUNCTION Real reduction_var<Var::beta>(REDUCE_FUNCTION_ARGS)
return ((gam - 1) * P(m_p.UU, k, j, i))/(0.5*(dot(Dtmp.bcon, Dtmp.bcov) + SMALL));
}

// Stuff that should be conserved
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::rhou0>(REDUCE_FUNCTION_ARGS)
{
return U(m_u.RHO, k, j, i);
}
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::mix_T00>(REDUCE_FUNCTION_ARGS)
{
return U(m_u.UU, k, j, i);
}
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::mix_T01>(REDUCE_FUNCTION_ARGS)
{
return U(m_u.U1, k, j, i);
}
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::mix_T02>(REDUCE_FUNCTION_ARGS)
{
return U(m_u.U2, k, j, i);
}
template <>
KOKKOS_INLINE_FUNCTION Real reduction_var<Var::mix_T03>(REDUCE_FUNCTION_ARGS)
{
return U(m_u.U3, k, j, i);
}

// Accretion rates: return a zone's contribution to the surface integral
// forming each rate measurement.
template <>
Expand Down
2 changes: 1 addition & 1 deletion pars/multizone/multizone.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ r_out = 134217728
r_in = 1

<parthenon/time>
tlim = 10000000000
tlim = 100000000000000000

<GRMHD>
cfl = 0.45
Expand Down

0 comments on commit 1b52648

Please sign in to comment.