Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

CalcMassFlux and CalculateMassAccretionRate ported to sparse packs #207

Closed
wants to merge 2 commits into from
Closed

CalcMassFlux and CalculateMassAccretionRate ported to sparse packs #207

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
500bc88
CalcMassFlux and CalculateMassAccretionRate ported to sparse packs
mari2895 Mar 14, 2024
9faa3ad
passing pack again
mari2895 Mar 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 14 additions & 17 deletions src/analysis/history.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,40 +26,37 @@ Real ReduceMassAccretionRate(MeshData<Real> *md) {
const auto jb = md->GetBoundsJ(IndexDomain::interior);
const auto kb = md->GetBoundsK(IndexDomain::interior);

namespace p = fluid_prim;
using parthenon::MakePackDescriptor;
Mesh *pmesh = md->GetMeshPointer();
auto &resolved_pkgs = pmesh->resolved_packages;
static auto desc = MakePackDescriptor<p::density, p::velocity>(resolved_pkgs.get());
auto v = desc.GetPack(md);
const int nblocks = v.GetNBlocks();

auto &pars = pmesh->packages.Get("geometry")->AllParams();
const Real xh = pars.Get<Real>("xh");

namespace p = fluid_prim;
const std::vector<std::string> vars({p::density::name(), p::velocity::name()});

PackIndexMap imap;
auto pack = md->PackVariables(vars, imap);

const int prho = imap[p::density::name()].first;
const int pvel_lo = imap[p::velocity::name()].first;
const int pvel_hi = imap[p::velocity::name()].second;

auto geom = Geometry::GetCoordinateSystem(md);

Real result = 0.0;
parthenon::par_reduce(
parthenon::LoopPatternMDRange(), "Phoebus History for Mass Accretion Rate",
DevExecSpace(), 0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
DevExecSpace(), 0, nblocks - 1, 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, Real &lresult) {
const auto &coords = pack.GetCoords(b);
const auto &coords = v.GetCoordinates(b);
if (coords.Xf<1>(i) <= xh && xh < coords.Xf<1>(i + 1)) {
const Real dx1 = coords.CellWidthFA(X1DIR, k, j, i);
const Real dx2 = coords.CellWidthFA(X2DIR, k, j, i);
const Real dx3 = coords.CellWidthFA(X3DIR, k, j, i);

// interp to make sure we're getting the horizon correct
auto m = (CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i + 1) -
CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i - 1)) /
auto m = (CalcMassFlux(md, v, geom, b, k, j, i + 1) -
CalcMassFlux(md, v, geom, b, k, j, i - 1)) /
(2.0 * dx1);
auto flux = (CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i) +
(xh - coords.Xc<1>(i)) * m) *
dx2 * dx3;
auto flux =
(CalcMassFlux(md, v, geom, b, k, j, i) + (xh - coords.Xc<1>(i)) * m) * dx2 *
dx3;

lresult += flux;
} else {
Expand Down
15 changes: 10 additions & 5 deletions src/analysis/history_utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "phoebus_utils/relativity_utils.hpp"
#include "phoebus_utils/robust.hpp"
#include "phoebus_utils/variables.hpp"
#include <interface/sparse_pack.hpp>
#include <kokkos_abstraction.hpp>
#include <parthenon/package.hpp>
#include <utils/error_checking.hpp>
Expand All @@ -32,25 +33,29 @@ using namespace Geometry;
namespace History {

template <typename Pack, typename Geometry>
KOKKOS_INLINE_FUNCTION Real CalcMassFlux(Pack &pack, Geometry &geom, const int prho,
const int pvel_lo, const int pvel_hi,
KOKKOS_INLINE_FUNCTION Real CalcMassFlux(MeshData<Real> *md, Pack &pack, Geometry &geom,
const int b, const int k, const int j,
const int i) {
namespace p = fluid_prim;
using parthenon::MakePackDescriptor;
Mesh *pmesh = md->GetMeshPointer();
auto &resolved_pkgs = pmesh->resolved_packages;
static auto desc = MakePackDescriptor<p::density, p::velocity>(resolved_pkgs.get());

Real gdet = geom.DetGamma(CellLocation::Cent, b, k, j, i);
Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i);
Real shift[3];
geom.ContravariantShift(CellLocation::Cent, b, k, j, i, shift);

const Real vel[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i),
pack(b, pvel_hi, k, j, i)};
const Real vel[] = {pack(b, p::velocity(), k, j, i), pack(b, p::velocity(1), k, j, i),
pack(b, p::velocity(2), k, j, i)};

Real gcov4[4][4];
geom.SpacetimeMetric(CellLocation::Cent, b, k, j, i, gcov4);
const Real W = phoebus::GetLorentzFactor(vel, gcov4);
const Real ucon = vel[0] - shift[0] * W / lapse;

return -lapse * gdet * pack(b, prho, k, j, i) * ucon;
return -lapse * gdet * pack(b, p::density(), k, j, i) * ucon;
}

template <typename Pack, typename Geometry>
Expand Down
Loading