From 4d0e6431f1025eb5f452ed2c2aa04877320edd69 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 4 Mar 2024 14:58:29 +0100 Subject: [PATCH 01/30] Add Phoebus tracers --- src/CMakeLists.txt | 2 + src/tracers/tracers.cpp | 294 ++++++++++++++++++++++++++++++++++++++++ src/tracers/tracers.hpp | 92 +++++++++++++ 3 files changed, 388 insertions(+) create mode 100644 src/tracers/tracers.cpp create mode 100644 src/tracers/tracers.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5d8fd375..c06b60ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,6 +18,8 @@ add_executable( hydro/srcterms/tabular_cooling.cpp refinement/gradient.cpp refinement/other.cpp + tracers/tracers.hpp + tracers/tracers.cpp utils/few_modes_ft.cpp ) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp new file mode 100644 index 00000000..99fd48dc --- /dev/null +++ b/src/tracers/tracers.cpp @@ -0,0 +1,294 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Tracer implementation refacored from https://github.com/lanl/phoebus +//======================================================================================== +// © 2021-2023. Triad National Security, LLC. All rights reserved. +// This program was produced under U.S. Government contract +// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +// is operated by Triad National Security, LLC for the U.S. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#include "tracers.hpp" +#include "geometry/geometry.hpp" +#include "geometry/geometry_utils.hpp" +#include "phoebus_utils/cell_locations.hpp" +#include "phoebus_utils/relativity_utils.hpp" +#include "phoebus_utils/variables.hpp" + +namespace tracers { +using namespace parthenon::package::prelude; + +std::shared_ptr Initialize(ParameterInput *pin) { + auto physics = std::make_shared("tracers"); + const bool active = pin->GetOrAddBoolean("physics", "tracers", false); + physics->AddParam("active", active); + if (!active) return physics; + + Params ¶ms = physics->AllParams(); + + const int num_tracers = pin->GetOrAddInteger("tracers", "num_tracers", 0); + params.Add("num_tracers", num_tracers); + + // Initialize random number generator pool + int rng_seed = pin->GetOrAddInteger("tracers", "rng_seed", time(NULL)); + physics->AddParam<>("rng_seed", rng_seed); + RNGPool rng_pool(rng_seed); + physics->AddParam<>("rng_pool", rng_pool); + + // Add swarm of tracers + std::string swarm_name = "tracers"; + Metadata swarm_metadata({Metadata::Provides, Metadata::None}); + physics->AddSwarm(swarm_name, swarm_metadata); + Metadata real_swarmvalue_metadata({Metadata::Real}); + physics->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); + + // thermo variables + physics->AddSwarmValue("rho", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("temperature", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("ye", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("entropy", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("pressure", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("energy", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("vel_x", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("vel_y", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("vel_z", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("lorentz", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("lapse", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("detgamma", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("shift_x", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("shift_y", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("shift_z", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("mass", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("bernoulli", swarm_name, real_swarmvalue_metadata); + + const bool mhd = pin->GetOrAddBoolean("fluid", "mhd", false); + + if (mhd) { + physics->AddSwarmValue("B_x", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("B_y", swarm_name, real_swarmvalue_metadata); + physics->AddSwarmValue("B_z", swarm_name, real_swarmvalue_metadata); + } + + return physics; +} // Initialize + +TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt) { + namespace p = fluid_prim; + + auto *pmb = rc->GetParentPointer(); + auto &sc = pmb->swarm_data.Get(); + auto &swarm = sc->Get("tracers"); + + const auto ndim = pmb->pmy_mesh->ndim; + + auto &x = swarm->Get("x").Get(); + auto &y = swarm->Get("y").Get(); + auto &z = swarm->Get("z").Get(); + + auto swarm_d = swarm->GetDeviceContext(); + + const std::vector vars = {p::velocity::name()}; + + PackIndexMap imap; + auto pack = rc->PackVariables(vars, imap); + + const int pvel_lo = imap[p::velocity::name()].first; + const int pvel_hi = imap[p::velocity::name()].second; + + auto geom = Geometry::GetCoordinateSystem(rc); + + // update loop. RK2 + const int max_active_index = swarm->GetMaxActiveIndex(); + pmb->par_for( + "Advect Tracers", 0, max_active_index, KOKKOS_LAMBDA(const int n) { + if (swarm_d.IsActive(n)) { + int k, j, i; + swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); + + Real rhs1, rhs2, rhs3; + + // predictor + tracers_rhs(pack, geom, pvel_lo, pvel_hi, ndim, dt, x(n), y(n), z(n), rhs1, + rhs2, rhs3); + const Real kx = x(n) + 0.5 * dt * rhs1; + const Real ky = y(n) + 0.5 * dt * rhs2; + const Real kz = z(n) + 0.5 * dt * rhs3; + + // corrector + tracers_rhs(pack, geom, pvel_lo, pvel_hi, ndim, dt, kx, ky, kz, rhs1, rhs2, + rhs3); + x(n) += rhs1 * dt; + y(n) += rhs2 * dt; + z(n) += rhs3 * dt; + + bool on_current_mesh_block = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + } + }); + + return TaskStatus::complete; +} // AdvectTracers + +/** + * FillDerived function for tracers. + * Registered Quantities (in addition to t, x, y, z): + * rho, T, ye, vel, energy, W_lorentz, pressure, + * lapse, shift, entropy, detgamma, B, bernoulli + **/ +void FillTracers(MeshBlockData *rc) { + using namespace LCInterp; + namespace p = fluid_prim; + + auto *pmb = rc->GetParentPointer(); + auto fluid = pmb->packages.Get("fluid"); + auto &sc = pmb->swarm_data.Get(); + auto &swarm = sc->Get("tracers"); + auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + + const auto mhd = fluid->Param("mhd"); + + // pull swarm vars + auto &x = swarm->Get("x").Get(); + auto &y = swarm->Get("y").Get(); + auto &z = swarm->Get("z").Get(); + auto &v1 = swarm->Get("vel_x").Get(); + auto &v2 = swarm->Get("vel_y").Get(); + auto &v3 = swarm->Get("vel_z").Get(); + auto B1 = v1.Get(); + auto B2 = v1.Get(); + auto B3 = v1.Get(); + if (mhd) { + B1 = swarm->Get("B_x").Get(); + B2 = swarm->Get("B_y").Get(); + B3 = swarm->Get("B_z").Get(); + } + auto &s_rho = swarm->Get("rho").Get(); + auto &s_temperature = swarm->Get("temperature").Get(); + auto &s_ye = swarm->Get("ye").Get(); + auto &s_entropy = swarm->Get("entropy").Get(); + auto &s_energy = swarm->Get("energy").Get(); + auto &s_lorentz = swarm->Get("lorentz").Get(); + auto &s_lapse = swarm->Get("lapse").Get(); + auto &s_shift_x = swarm->Get("shift_x").Get(); + auto &s_shift_y = swarm->Get("shift_y").Get(); + auto &s_shift_z = swarm->Get("shift_z").Get(); + auto &s_detgamma = swarm->Get("detgamma").Get(); + auto &s_pressure = swarm->Get("pressure").Get(); + auto &s_bernoulli = swarm->Get("bernoulli").Get(); + + auto swarm_d = swarm->GetDeviceContext(); + + std::vector vars = {p::density::name(), p::temperature::name(), + p::velocity::name(), p::energy::name(), + p::pressure::name()}; + if (mhd) { + vars.push_back(p::bfield::name()); + } + + PackIndexMap imap; + auto pack = rc->PackVariables(vars, imap); + + const int pvel_lo = imap[p::velocity::name()].first; + const int pvel_hi = imap[p::velocity::name()].second; + const int pB_lo = imap[p::bfield::name()].first; + const int pB_hi = imap[p::bfield::name()].second; + const int prho = imap[p::density::name()].first; + const int ptemp = imap[p::temperature::name()].first; + const int pye = imap[p::ye::name()].second; + const int penergy = imap[p::energy::name()].first; + const int ppres = imap[p::pressure::name()].first; + + auto geom = Geometry::GetCoordinateSystem(rc); + // update loop. + const int max_active_index = swarm->GetMaxActiveIndex(); + pmb->par_for( + "Fill Tracers", 0, max_active_index, KOKKOS_LAMBDA(const int n) { + if (swarm_d.IsActive(n)) { + int k, j, i; + swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); + + // geom quantities + Real gcov4[4][4]; + geom.SpacetimeMetric(0.0, x(n), y(n), z(n), gcov4); + Real lapse = geom.Lapse(0.0, x(n), y(n), z(n)); + Real shift[3]; + geom.ContravariantShift(0.0, x(n), y(n), z(n), shift); + const Real gdet = geom.DetGamma(0.0, x(n), y(n), z(n)); + + // Interpolate + const Real Wvel_X1 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_lo); + const Real Wvel_X2 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_lo + 1); + const Real Wvel_X3 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_hi); + Real B_X1 = 0.0; + Real B_X2 = 0.0; + Real B_X3 = 0.0; + if (mhd) { + B_X1 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_lo); + B_X2 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_lo + 1); + B_X3 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_hi); + } + const Real rho = LCInterp::Do(0, x(n), y(n), z(n), pack, prho); + const Real temperature = LCInterp::Do(0, x(n), y(n), z(n), pack, ptemp); + const Real energy = LCInterp::Do(0, x(n), y(n), z(n), pack, penergy); + const Real pressure = LCInterp::Do(0, x(n), y(n), z(n), pack, ppres); + const Real Wvel[] = {Wvel_X1, Wvel_X2, Wvel_X3}; + const Real W = phoebus::GetLorentzFactor(Wvel, gcov4); + const Real vel_X1 = Wvel_X1 / W; + const Real vel_X2 = Wvel_X2 / W; + const Real vel_X3 = Wvel_X3 / W; + Real ye; + Real lambda[2] = {0.0, 0.0}; + if (pye > 0) { + ye = LCInterp::Do(0, x(n), y(n), z(n), pack, pye); + lambda[1] = ye; + } else { + ye = 0.0; + } + const Real entropy = + eos.EntropyFromDensityTemperature(rho, temperature, lambda); + + // bernoulli + const Real h = 1.0 + energy + pressure / rho; + const Real bernoulli = -(W / lapse) * h - 1.0; + + // store + s_rho(n) = rho; + s_temperature(n) = temperature; + s_ye(n) = ye; + s_energy(n) = energy; + s_entropy(n) = entropy; + v1(n) = vel_X1; + v2(n) = vel_X3; + v3(n) = vel_X2; + s_shift_x(n) = shift[0]; + s_shift_y(n) = shift[1]; + s_shift_z(n) = shift[2]; + s_lapse(n) = lapse; + s_lorentz(n) = W; + s_detgamma(n) = gdet; + s_pressure(n) = pressure; + s_bernoulli(n) = bernoulli; + if (mhd) { + B1(n) = B_X1; + B2(n) = B_X2; + B3(n) = B_X3; + } + + bool on_current_mesh_block = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + } + }); + +} // FillTracers + +} // namespace tracers diff --git a/src/tracers/tracers.hpp b/src/tracers/tracers.hpp new file mode 100644 index 00000000..9ebbd962 --- /dev/null +++ b/src/tracers/tracers.hpp @@ -0,0 +1,92 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Tracer implementation refacored from https://github.com/lanl/phoebus +//======================================================================================== +// © 2021-2023. Triad National Security, LLC. All rights reserved. +// This program was produced under U.S. Government contract +// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +// is operated by Triad National Security, LLC for the U.S. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#ifndef TRACERS_HPP_ +#define TRACERS_HPP_ + +#include + +#include "Kokkos_Random.hpp" + +#include +#include + +#include "geometry/geometry.hpp" +#include "geometry/geometry_utils.hpp" +#include "microphysics/eos_phoebus/eos_phoebus.hpp" +#include "phoebus_utils/cell_locations.hpp" +#include "phoebus_utils/phoebus_interpolation.hpp" +#include "phoebus_utils/relativity_utils.hpp" +#include "phoebus_utils/variables.hpp" + +using namespace parthenon::driver::prelude; +using namespace parthenon::package::prelude; +using namespace parthenon; +using Microphysics::EOS::EOS; + +typedef Kokkos::Random_XorShift64_Pool<> RNGPool; + +namespace tracers { + +std::shared_ptr Initialize(ParameterInput *pin); + +/** + * RHS of tracer advection equations + * alpha v^i - beta^i + * dt not included + **/ +template +KOKKOS_INLINE_FUNCTION void tracers_rhs(Pack &pack, Geometry &geom, const int pvel_lo, + const int pvel_hi, const int ndim, const Real dt, + const Real x, const Real y, const Real z, + Real &rhs1, Real &rhs2, Real &rhs3) { + + // geom quantities + Real gcov4[4][4]; + geom.SpacetimeMetric(0.0, x, y, z, gcov4); + Real lapse = geom.Lapse(0.0, x, y, z); + Real shift[3]; + geom.ContravariantShift(0.0, x, y, z, shift); + + // Get shift, W, lapse + const Real Wvel_X1 = LCInterp::Do(0, x, y, z, pack, pvel_lo); + const Real Wvel_X2 = LCInterp::Do(0, x, y, z, pack, pvel_lo + 1); + const Real Wvel_X3 = LCInterp::Do(0, x, y, z, pack, pvel_hi); + const Real Wvel[] = {Wvel_X1, Wvel_X2, Wvel_X3}; + const Real W = phoebus::GetLorentzFactor(Wvel, gcov4); + const Real vel_X1 = Wvel_X1 / W; + const Real vel_X2 = Wvel_X2 / W; + const Real vel_X3 = Wvel_X3 / W; + + rhs1 = (lapse * vel_X1 - shift[0]); + rhs2 = (lapse * vel_X2 - shift[1]); + rhs3 = 0.0; + if (ndim > 2) { + rhs3 = (lapse * vel_X3 - shift[2]); + } +} + +TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt); + +void FillTracers(MeshBlockData *rc); + +} // namespace tracers + +#endif // TRACERS_HPP_ From 4c5a77760f242d9b2226cb7b4a95091dc2972893 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 4 Mar 2024 20:00:09 +0100 Subject: [PATCH 02/30] Fix/adapt tracers to AthenaPK --- src/CMakeLists.txt | 2 +- src/hydro/hydro.cpp | 2 + src/hydro/hydro_driver.cpp | 31 +++++ src/tracers/tracers.cpp | 270 +++++++++++-------------------------- src/tracers/tracers.hpp | 54 +------- 5 files changed, 119 insertions(+), 240 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c06b60ec..40c5eb40 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,8 +18,8 @@ add_executable( hydro/srcterms/tabular_cooling.cpp refinement/gradient.cpp refinement/other.cpp - tracers/tracers.hpp tracers/tracers.cpp + tracers/tracers.hpp utils/few_modes_ft.cpp ) diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index e5402454..6dff2e72 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -24,6 +24,7 @@ #include "../recon/weno3_simple.hpp" #include "../recon/wenoz_simple.hpp" #include "../refinement/refinement.hpp" +#include "../tracers/tracers.hpp" #include "../units.hpp" #include "defs.hpp" #include "diffusion/diffusion.hpp" @@ -52,6 +53,7 @@ using parthenon::HistoryOutputVar; parthenon::Packages_t ProcessPackages(std::unique_ptr &pin) { parthenon::Packages_t packages; packages.Add(Hydro::Initialize(pin.get())); + packages.Add(Tracers::Initialize(pin.get())); return packages; } diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index e33b1b0d..84e1778d 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -19,6 +19,7 @@ #include "../eos/adiabatic_hydro.hpp" #include "../pgen/cluster/agn_triggering.hpp" #include "../pgen/cluster/magnetic_tower.hpp" +#include "../tracers/tracers.hpp" #include "glmmhd/glmmhd.hpp" #include "hydro.hpp" #include "hydro_driver.hpp" @@ -382,6 +383,36 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { fill_derived, parthenon::Update::EstimateTimestep>, mu0.get()); } } + auto tracers_pkg = pmesh->packages.Get("tracers"); + // First order operator split tracer advection + if (stage == integrator->nstages && tracers_pkg->Param("enabled")) { + const std::string swarm_name = "tracers"; + TaskRegion &sync_region_tr = tc.AddRegion(1); + { + for (auto &pmb : blocks) { + auto &tl = sync_region_tr[0]; + auto &sd = pmb->swarm_data.Get(); + auto reset_comms = + tl.AddTask(none, &SwarmContainer::ResetCommunication, sd.get()); + } + } + + TaskRegion &async_region_tr = tc.AddRegion(blocks.size()); + for (int n = 0; n < blocks.size(); n++) { + auto &tl = async_region_tr[n]; + auto &pmb = blocks[n]; + auto &sd = pmb->swarm_data.Get(); + auto &mbd0 = pmb->meshblock_data.Get("base"); + auto tracer_advect = + tl.AddTask(none, Tracers::AdvectTracers, mbd0.get(), integrator->dt); + + auto send = tl.AddTask(tracer_advect, &SwarmContainer::Send, sd.get(), + BoundaryCommSubset::all); + + auto receive = + tl.AddTask(send, &SwarmContainer::Receive, sd.get(), BoundaryCommSubset::all); + } + } if (stage == integrator->nstages && pmesh->adaptive) { TaskRegion &async_region_4 = tc.AddRegion(num_task_lists_executed_independently); diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 99fd48dc..2f5bad1a 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -19,75 +19,65 @@ // publicly, and to permit others to do so. #include "tracers.hpp" -#include "geometry/geometry.hpp" -#include "geometry/geometry_utils.hpp" -#include "phoebus_utils/cell_locations.hpp" -#include "phoebus_utils/relativity_utils.hpp" -#include "phoebus_utils/variables.hpp" +#include "../main.hpp" +#include "utils/error_checking.hpp" -namespace tracers { +namespace Tracers { using namespace parthenon::package::prelude; std::shared_ptr Initialize(ParameterInput *pin) { - auto physics = std::make_shared("tracers"); - const bool active = pin->GetOrAddBoolean("physics", "tracers", false); - physics->AddParam("active", active); - if (!active) return physics; + auto tracer_pkg = std::make_shared("tracers"); + const bool enabled = pin->GetOrAddBoolean("tracers", "enabled", false); + tracer_pkg->AddParam("enabled", enabled); + if (!enabled) return tracer_pkg; - Params ¶ms = physics->AllParams(); + Params ¶ms = tracer_pkg->AllParams(); const int num_tracers = pin->GetOrAddInteger("tracers", "num_tracers", 0); params.Add("num_tracers", num_tracers); // Initialize random number generator pool int rng_seed = pin->GetOrAddInteger("tracers", "rng_seed", time(NULL)); - physics->AddParam<>("rng_seed", rng_seed); + tracer_pkg->AddParam<>("rng_seed", rng_seed); RNGPool rng_pool(rng_seed); - physics->AddParam<>("rng_pool", rng_pool); + tracer_pkg->AddParam<>("rng_pool", rng_pool); // Add swarm of tracers std::string swarm_name = "tracers"; Metadata swarm_metadata({Metadata::Provides, Metadata::None}); - physics->AddSwarm(swarm_name, swarm_metadata); + tracer_pkg->AddSwarm(swarm_name, swarm_metadata); Metadata real_swarmvalue_metadata({Metadata::Real}); - physics->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); + tracer_pkg->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); + // TODO(pgrete) Add CheckDesired/required for vars // thermo variables - physics->AddSwarmValue("rho", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("temperature", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("ye", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("entropy", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("pressure", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("energy", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("vel_x", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("vel_y", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("vel_z", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("lorentz", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("lapse", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("detgamma", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("shift_x", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("shift_y", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("shift_z", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("mass", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("bernoulli", swarm_name, real_swarmvalue_metadata); - - const bool mhd = pin->GetOrAddBoolean("fluid", "mhd", false); + tracer_pkg->AddSwarmValue("rho", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("pressure", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("vel_x", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("vel_y", swarm_name, real_swarmvalue_metadata); + // TODO(pgrete) check proper handling of <3D sims + tracer_pkg->AddSwarmValue("vel_z", swarm_name, real_swarmvalue_metadata); + + // TODO(pgrete) this should be safe because we call this package init after the hydro + // one, but we should check if there's direct way to access Params of other packages. + const bool mhd = pin->GetString("hydro", "fluid") == "glmmhd"; + + PARTHENON_REQUIRE_THROWS(pin->GetString("parthenon/mesh", "refinement") == "none", + "Tracers/swarms currently only supported on uniform meshes."); if (mhd) { - physics->AddSwarmValue("B_x", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("B_y", swarm_name, real_swarmvalue_metadata); - physics->AddSwarmValue("B_z", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("B_x", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("B_y", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("B_z", swarm_name, real_swarmvalue_metadata); } - return physics; + return tracer_pkg; } // Initialize -TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt) { - namespace p = fluid_prim; - - auto *pmb = rc->GetParentPointer(); - auto &sc = pmb->swarm_data.Get(); - auto &swarm = sc->Get("tracers"); +TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { + auto *pmb = mbd->GetParentPointer(); + auto &sd = pmb->swarm_data.Get(); + auto &swarm = sd->Get("tracers"); const auto ndim = pmb->pmy_mesh->ndim; @@ -97,15 +87,7 @@ TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt) { auto swarm_d = swarm->GetDeviceContext(); - const std::vector vars = {p::velocity::name()}; - - PackIndexMap imap; - auto pack = rc->PackVariables(vars, imap); - - const int pvel_lo = imap[p::velocity::name()].first; - const int pvel_hi = imap[p::velocity::name()].second; - - auto geom = Geometry::GetCoordinateSystem(rc); + const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); // update loop. RK2 const int max_active_index = swarm->GetMaxActiveIndex(); @@ -115,24 +97,20 @@ TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt) { int k, j, i; swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); - Real rhs1, rhs2, rhs3; - + // Predictor/corrector will first make sense with non constant interpolation + // TODO(pgrete) add non-cnonst interpolation // predictor - tracers_rhs(pack, geom, pvel_lo, pvel_hi, ndim, dt, x(n), y(n), z(n), rhs1, - rhs2, rhs3); - const Real kx = x(n) + 0.5 * dt * rhs1; - const Real ky = y(n) + 0.5 * dt * rhs2; - const Real kz = z(n) + 0.5 * dt * rhs3; + // const Real kx = x(n) + 0.5 * dt * rhs1; + // const Real ky = y(n) + 0.5 * dt * rhs2; + // const Real kz = z(n) + 0.5 * dt * rhs3; // corrector - tracers_rhs(pack, geom, pvel_lo, pvel_hi, ndim, dt, kx, ky, kz, rhs1, rhs2, - rhs3); - x(n) += rhs1 * dt; - y(n) += rhs2 * dt; - z(n) += rhs3 * dt; - - bool on_current_mesh_block = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + x(n) += prim_pack(IV1, k, j, i) * dt; + y(n) += prim_pack(IV2, k, j, i) * dt; + z(n) += prim_pack(IV3, k, j, i) * dt; + + bool unused_temp = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unused_temp); } }); @@ -142,73 +120,42 @@ TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt) { /** * FillDerived function for tracers. * Registered Quantities (in addition to t, x, y, z): - * rho, T, ye, vel, energy, W_lorentz, pressure, - * lapse, shift, entropy, detgamma, B, bernoulli + * rho, vel, B **/ -void FillTracers(MeshBlockData *rc) { - using namespace LCInterp; - namespace p = fluid_prim; +void FillTracers(MeshBlockData *mbd) { - auto *pmb = rc->GetParentPointer(); - auto fluid = pmb->packages.Get("fluid"); - auto &sc = pmb->swarm_data.Get(); - auto &swarm = sc->Get("tracers"); - auto eos = pmb->packages.Get("eos")->Param("d.EOS"); + auto *pmb = mbd->GetParentPointer(); + auto hydro_pkg = pmb->packages.Get("Hydro"); + auto &sd = pmb->swarm_data.Get(); + auto &swarm = sd->Get("tracers"); - const auto mhd = fluid->Param("mhd"); + const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; + // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) // pull swarm vars auto &x = swarm->Get("x").Get(); auto &y = swarm->Get("y").Get(); auto &z = swarm->Get("z").Get(); - auto &v1 = swarm->Get("vel_x").Get(); - auto &v2 = swarm->Get("vel_y").Get(); - auto &v3 = swarm->Get("vel_z").Get(); - auto B1 = v1.Get(); - auto B2 = v1.Get(); - auto B3 = v1.Get(); + auto &vel_x = swarm->Get("vel_x").Get(); + auto &vel_y = swarm->Get("vel_y").Get(); + auto &vel_z = swarm->Get("vel_z").Get(); + // Assign some (definitely existing) default var + auto B_x = vel_x.Get(); + auto B_y = vel_x.Get(); + auto B_z = vel_x.Get(); if (mhd) { - B1 = swarm->Get("B_x").Get(); - B2 = swarm->Get("B_y").Get(); - B3 = swarm->Get("B_z").Get(); + B_x = swarm->Get("B_x").Get(); + B_y = swarm->Get("B_y").Get(); + B_z = swarm->Get("B_z").Get(); } - auto &s_rho = swarm->Get("rho").Get(); - auto &s_temperature = swarm->Get("temperature").Get(); - auto &s_ye = swarm->Get("ye").Get(); - auto &s_entropy = swarm->Get("entropy").Get(); - auto &s_energy = swarm->Get("energy").Get(); - auto &s_lorentz = swarm->Get("lorentz").Get(); - auto &s_lapse = swarm->Get("lapse").Get(); - auto &s_shift_x = swarm->Get("shift_x").Get(); - auto &s_shift_y = swarm->Get("shift_y").Get(); - auto &s_shift_z = swarm->Get("shift_z").Get(); - auto &s_detgamma = swarm->Get("detgamma").Get(); - auto &s_pressure = swarm->Get("pressure").Get(); - auto &s_bernoulli = swarm->Get("bernoulli").Get(); + auto &rho = swarm->Get("rho").Get(); + auto &pressure = swarm->Get("pressure").Get(); - auto swarm_d = swarm->GetDeviceContext(); + // Get hydro/mhd fluid vars + const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); - std::vector vars = {p::density::name(), p::temperature::name(), - p::velocity::name(), p::energy::name(), - p::pressure::name()}; - if (mhd) { - vars.push_back(p::bfield::name()); - } - - PackIndexMap imap; - auto pack = rc->PackVariables(vars, imap); - - const int pvel_lo = imap[p::velocity::name()].first; - const int pvel_hi = imap[p::velocity::name()].second; - const int pB_lo = imap[p::bfield::name()].first; - const int pB_hi = imap[p::bfield::name()].second; - const int prho = imap[p::density::name()].first; - const int ptemp = imap[p::temperature::name()].first; - const int pye = imap[p::ye::name()].second; - const int penergy = imap[p::energy::name()].first; - const int ppres = imap[p::pressure::name()].first; + auto swarm_d = swarm->GetDeviceContext(); - auto geom = Geometry::GetCoordinateSystem(rc); // update loop. const int max_active_index = swarm->GetMaxActiveIndex(); pmb->par_for( @@ -217,78 +164,23 @@ void FillTracers(MeshBlockData *rc) { int k, j, i; swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); - // geom quantities - Real gcov4[4][4]; - geom.SpacetimeMetric(0.0, x(n), y(n), z(n), gcov4); - Real lapse = geom.Lapse(0.0, x(n), y(n), z(n)); - Real shift[3]; - geom.ContravariantShift(0.0, x(n), y(n), z(n), shift); - const Real gdet = geom.DetGamma(0.0, x(n), y(n), z(n)); - - // Interpolate - const Real Wvel_X1 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_lo); - const Real Wvel_X2 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_lo + 1); - const Real Wvel_X3 = LCInterp::Do(0, x(n), y(n), z(n), pack, pvel_hi); - Real B_X1 = 0.0; - Real B_X2 = 0.0; - Real B_X3 = 0.0; - if (mhd) { - B_X1 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_lo); - B_X2 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_lo + 1); - B_X3 = LCInterp::Do(0, x(n), y(n), z(n), pack, pB_hi); - } - const Real rho = LCInterp::Do(0, x(n), y(n), z(n), pack, prho); - const Real temperature = LCInterp::Do(0, x(n), y(n), z(n), pack, ptemp); - const Real energy = LCInterp::Do(0, x(n), y(n), z(n), pack, penergy); - const Real pressure = LCInterp::Do(0, x(n), y(n), z(n), pack, ppres); - const Real Wvel[] = {Wvel_X1, Wvel_X2, Wvel_X3}; - const Real W = phoebus::GetLorentzFactor(Wvel, gcov4); - const Real vel_X1 = Wvel_X1 / W; - const Real vel_X2 = Wvel_X2 / W; - const Real vel_X3 = Wvel_X3 / W; - Real ye; - Real lambda[2] = {0.0, 0.0}; - if (pye > 0) { - ye = LCInterp::Do(0, x(n), y(n), z(n), pack, pye); - lambda[1] = ye; - } else { - ye = 0.0; - } - const Real entropy = - eos.EntropyFromDensityTemperature(rho, temperature, lambda); - - // bernoulli - const Real h = 1.0 + energy + pressure / rho; - const Real bernoulli = -(W / lapse) * h - 1.0; - - // store - s_rho(n) = rho; - s_temperature(n) = temperature; - s_ye(n) = ye; - s_energy(n) = energy; - s_entropy(n) = entropy; - v1(n) = vel_X1; - v2(n) = vel_X3; - v3(n) = vel_X2; - s_shift_x(n) = shift[0]; - s_shift_y(n) = shift[1]; - s_shift_z(n) = shift[2]; - s_lapse(n) = lapse; - s_lorentz(n) = W; - s_detgamma(n) = gdet; - s_pressure(n) = pressure; - s_bernoulli(n) = bernoulli; + // TODO(pgrete) Interpolate + rho(n) = prim_pack(IDN, k, j, i); + vel_x(n) = prim_pack(IV1, k, j, i); + vel_y(n) = prim_pack(IV2, k, j, i); + vel_z(n) = prim_pack(IV3, k, j, i); + pressure(n) = prim_pack(IPR, k, j, i); if (mhd) { - B1(n) = B_X1; - B2(n) = B_X2; - B3(n) = B_X3; + B_x(n) = prim_pack(IB1, k, j, i); + B_y(n) = prim_pack(IB2, k, j, i); + B_z(n) = prim_pack(IB3, k, j, i); } - bool on_current_mesh_block = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + bool unsed_tmp = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } }); } // FillTracers -} // namespace tracers +} // namespace Tracers diff --git a/src/tracers/tracers.hpp b/src/tracers/tracers.hpp index 9ebbd962..1a21b2a7 100644 --- a/src/tracers/tracers.hpp +++ b/src/tracers/tracers.hpp @@ -28,65 +28,19 @@ #include #include -#include "geometry/geometry.hpp" -#include "geometry/geometry_utils.hpp" -#include "microphysics/eos_phoebus/eos_phoebus.hpp" -#include "phoebus_utils/cell_locations.hpp" -#include "phoebus_utils/phoebus_interpolation.hpp" -#include "phoebus_utils/relativity_utils.hpp" -#include "phoebus_utils/variables.hpp" - using namespace parthenon::driver::prelude; using namespace parthenon::package::prelude; -using namespace parthenon; -using Microphysics::EOS::EOS; -typedef Kokkos::Random_XorShift64_Pool<> RNGPool; +using RNGPool = Kokkos::Random_XorShift64_Pool<>; -namespace tracers { +namespace Tracers { std::shared_ptr Initialize(ParameterInput *pin); -/** - * RHS of tracer advection equations - * alpha v^i - beta^i - * dt not included - **/ -template -KOKKOS_INLINE_FUNCTION void tracers_rhs(Pack &pack, Geometry &geom, const int pvel_lo, - const int pvel_hi, const int ndim, const Real dt, - const Real x, const Real y, const Real z, - Real &rhs1, Real &rhs2, Real &rhs3) { - - // geom quantities - Real gcov4[4][4]; - geom.SpacetimeMetric(0.0, x, y, z, gcov4); - Real lapse = geom.Lapse(0.0, x, y, z); - Real shift[3]; - geom.ContravariantShift(0.0, x, y, z, shift); - - // Get shift, W, lapse - const Real Wvel_X1 = LCInterp::Do(0, x, y, z, pack, pvel_lo); - const Real Wvel_X2 = LCInterp::Do(0, x, y, z, pack, pvel_lo + 1); - const Real Wvel_X3 = LCInterp::Do(0, x, y, z, pack, pvel_hi); - const Real Wvel[] = {Wvel_X1, Wvel_X2, Wvel_X3}; - const Real W = phoebus::GetLorentzFactor(Wvel, gcov4); - const Real vel_X1 = Wvel_X1 / W; - const Real vel_X2 = Wvel_X2 / W; - const Real vel_X3 = Wvel_X3 / W; - - rhs1 = (lapse * vel_X1 - shift[0]); - rhs2 = (lapse * vel_X2 - shift[1]); - rhs3 = 0.0; - if (ndim > 2) { - rhs3 = (lapse * vel_X3 - shift[2]); - } -} - -TaskStatus AdvectTracers(MeshBlockData *rc, const Real dt); +TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt); void FillTracers(MeshBlockData *rc); -} // namespace tracers +} // namespace Tracers #endif // TRACERS_HPP_ From b993570c4909be509ba4919abaa9bc8758b68a4a Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 4 Mar 2024 20:00:20 +0100 Subject: [PATCH 03/30] Add tracer support to turbulence pgen --- src/pgen/turbulence.cpp | 61 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 2c9ff9a7..bd600261 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -26,6 +26,7 @@ // AthenaPK headers #include "../main.hpp" +#include "../tracers/tracers.hpp" #include "../units.hpp" #include "../utils/few_modes_ft.hpp" @@ -333,6 +334,59 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) + SQR(u(IB3, k, j, i))); } }); + + /* tracer init section */ + auto tracer_pkg = pmb->packages.Get("tracers"); + if (tracer_pkg->Param("enabled")) { + + for (int b = 0; b < md->NumBlocks(); b++) { + pmb = md->GetBlockData(0)->GetBlockPointer(); + auto &sd = pmb->swarm_data.Get(); + auto &swarm = pmb->swarm_data.Get()->Get("tracers"); + auto rng_pool = tracer_pkg->Param("rng_pool"); + + const Real &x_min = pmb->coords.Xf<1>(ib.s); + const Real &y_min = pmb->coords.Xf<2>(jb.s); + const Real &z_min = pmb->coords.Xf<3>(kb.s); + const Real &x_max = pmb->coords.Xf<1>(ib.e + 1); + const Real &y_max = pmb->coords.Xf<2>(jb.e + 1); + const Real &z_max = pmb->coords.Xf<3>(kb.e + 1); + + // as for num_tracers on each block... will get too many on multiple blocks + // TODO: distribute amongst blocks. + const auto num_tracers_total = tracer_pkg->Param("num_tracers"); + const int number_block = num_tracers_total; + + ParArrayND new_indices; + swarm->AddEmptyParticles(number_block, new_indices); + + auto &x = swarm->Get("x").Get(); + auto &y = swarm->Get("y").Get(); + auto &z = swarm->Get("z").Get(); + auto &id = swarm->Get("id").Get(); + + auto swarm_d = swarm->GetDeviceContext(); + + const int gid = pmb->gid; + const int max_active_index = swarm->GetMaxActiveIndex(); + pmb->par_for( + "ProblemGenerator::Turbulence::DistributeTracers", 0, max_active_index, + KOKKOS_LAMBDA(const int n) { + if (swarm_d.IsActive(n)) { + auto rng_gen = rng_pool.get_state(); + + x(n) = x_min + rng_gen.drand() * (x_max - x_min); + y(n) = y_min + rng_gen.drand() * (y_max - y_min); + z(n) = z_min + rng_gen.drand() * (z_max - z_min); + id(n) = num_tracers_total * gid + n; + + bool on_current_mesh_block = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); + rng_pool.free_state(rng_gen); + } + }); + } + } } //---------------------------------------------------------------------------------------- @@ -476,6 +530,13 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { // store state of distribution auto state_dist = few_modes_ft.GetDistState(); pin->SetString("problem/turbulence", "state_dist", state_dist); + + // TODO(pgrete) move this to a more central location in the driver + auto tracer_pkg = pmb->packages.Get("tracers"); + if (tracer_pkg->Param("enabled")) { + auto &mbd = pmb->meshblock_data.Get(); + Tracers::FillTracers(mbd.get()); + } } } // namespace turbulence From 1465cd487830157cc1185b97df75f218f1930558 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 6 Mar 2024 11:54:35 +0100 Subject: [PATCH 04/30] Fix swarm interface --- external/parthenon | 2 +- src/pgen/turbulence.cpp | 43 +++++++++++++++++++++++------------------ src/tracers/tracers.cpp | 14 +++++++++----- 3 files changed, 34 insertions(+), 25 deletions(-) diff --git a/external/parthenon b/external/parthenon index 88ece809..3225612e 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 88ece809dc56ae16d695db8d9051433bff01feb1 +Subproject commit 3225612e267456e5731f80c594d45e6457274c3f diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index bd600261..484f69b3 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -354,11 +354,12 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { // as for num_tracers on each block... will get too many on multiple blocks // TODO: distribute amongst blocks. - const auto num_tracers_total = tracer_pkg->Param("num_tracers"); - const int number_block = num_tracers_total; + const auto num_tracers_per_cell = tracer_pkg->Param("num_tracers_per_cell"); + const auto num_tracers_per_block = + static_cast(pmb->GetNumberOfMeshBlockCells() * num_tracers_per_cell); - ParArrayND new_indices; - swarm->AddEmptyParticles(number_block, new_indices); + // Create new particles and get accessor + auto new_particles_context = swarm->AddEmptyParticles(num_tracers_per_block); auto &x = swarm->Get("x").Get(); auto &y = swarm->Get("y").Get(); @@ -368,22 +369,26 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { auto swarm_d = swarm->GetDeviceContext(); const int gid = pmb->gid; - const int max_active_index = swarm->GetMaxActiveIndex(); pmb->par_for( - "ProblemGenerator::Turbulence::DistributeTracers", 0, max_active_index, - KOKKOS_LAMBDA(const int n) { - if (swarm_d.IsActive(n)) { - auto rng_gen = rng_pool.get_state(); - - x(n) = x_min + rng_gen.drand() * (x_max - x_min); - y(n) = y_min + rng_gen.drand() * (y_max - y_min); - z(n) = z_min + rng_gen.drand() * (z_max - z_min); - id(n) = num_tracers_total * gid + n; - - bool on_current_mesh_block = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); - rng_pool.free_state(rng_gen); - } + "ProblemGenerator::Turbulence::DistributeTracers", 0, + new_particles_context.GetNewParticlesMaxIndex(), + KOKKOS_LAMBDA(const int new_n) { + auto rng_gen = rng_pool.get_state(); + const int n = new_particles_context.GetNewParticleIndex(new_n); + + x(n) = x_min + rng_gen.drand() * (x_max - x_min); + y(n) = y_min + rng_gen.drand() * (y_max - y_min); + z(n) = z_min + rng_gen.drand() * (z_max - z_min); + // Note that his works only during one time init. + // If (somehwere else) we eventually add dynamic particles, then we need to + // manage ids (not indices) more globally. + id(n) = num_tracers_per_block * gid + n; + + rng_pool.free_state(rng_gen); + + // TODO(pgrete) check if this actually required + bool on_current_mesh_block = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block); }); } } diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 2f5bad1a..0cf5f204 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -33,8 +33,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { Params ¶ms = tracer_pkg->AllParams(); - const int num_tracers = pin->GetOrAddInteger("tracers", "num_tracers", 0); - params.Add("num_tracers", num_tracers); + const auto num_tracers_per_cell = + pin->GetOrAddReal("tracers", "num_tracers_per_cell", 0.0); + PARTHENON_REQUIRE_THROWS(num_tracers_per_cell >= 0.0, + "What's a negative number of particles per cell?!"); + params.Add("num_tracers_per_cell", num_tracers_per_cell); // Initialize random number generator pool int rng_seed = pin->GetOrAddInteger("tracers", "rng_seed", time(NULL)); @@ -57,12 +60,15 @@ std::shared_ptr Initialize(ParameterInput *pin) { tracer_pkg->AddSwarmValue("vel_y", swarm_name, real_swarmvalue_metadata); // TODO(pgrete) check proper handling of <3D sims tracer_pkg->AddSwarmValue("vel_z", swarm_name, real_swarmvalue_metadata); + PARTHENON_REQUIRE_THROWS(pin->GetInteger("parthenon/mesh", "nx3") > 1, + "Tracers/swarms currently only supported/tested in 3D."); // TODO(pgrete) this should be safe because we call this package init after the hydro // one, but we should check if there's direct way to access Params of other packages. const bool mhd = pin->GetString("hydro", "fluid") == "glmmhd"; - PARTHENON_REQUIRE_THROWS(pin->GetString("parthenon/mesh", "refinement") == "none", + PARTHENON_REQUIRE_THROWS(pin->GetOrAddString("parthenon/mesh", "refinement", "none") == + "none", "Tracers/swarms currently only supported on uniform meshes."); if (mhd) { @@ -79,8 +85,6 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { auto &sd = pmb->swarm_data.Get(); auto &swarm = sd->Get("tracers"); - const auto ndim = pmb->pmy_mesh->ndim; - auto &x = swarm->Get("x").Get(); auto &y = swarm->Get("y").Get(); auto &z = swarm->Get("z").Get(); From db3c876fd61e0f68b96fe519cb4911bb4b76f265 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 6 Mar 2024 15:21:46 +0100 Subject: [PATCH 05/30] Fix particle init --- src/pgen/turbulence.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 484f69b3..e2c833d5 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -340,7 +340,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { if (tracer_pkg->Param("enabled")) { for (int b = 0; b < md->NumBlocks(); b++) { - pmb = md->GetBlockData(0)->GetBlockPointer(); + pmb = md->GetBlockData(b)->GetBlockPointer(); auto &sd = pmb->swarm_data.Get(); auto &swarm = pmb->swarm_data.Get()->Get("tracers"); auto rng_pool = tracer_pkg->Param("rng_pool"); From 13505df722b7b745a1625cf499d1e4633011ab8b Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 6 Mar 2024 21:58:17 +0100 Subject: [PATCH 06/30] Call FillTracer every cycle --- src/hydro/hydro_driver.cpp | 3 +++ src/pgen/turbulence.cpp | 7 ------- src/tracers/tracers.cpp | 4 +++- src/tracers/tracers.hpp | 2 +- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index 84e1778d..b352c382 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -12,6 +12,7 @@ // Parthenon headers #include "amr_criteria/refinement_package.hpp" +#include "basic_types.hpp" #include "bvals/comms/bvals_in_one.hpp" #include "prolong_restrict/prolong_restrict.hpp" #include @@ -411,6 +412,8 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto receive = tl.AddTask(send, &SwarmContainer::Receive, sd.get(), BoundaryCommSubset::all); + + auto fill = tl.AddTask(receive, Tracers::FillTracers, mbd0.get(), tm); } } diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index e2c833d5..23f101ee 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -535,13 +535,6 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { // store state of distribution auto state_dist = few_modes_ft.GetDistState(); pin->SetString("problem/turbulence", "state_dist", state_dist); - - // TODO(pgrete) move this to a more central location in the driver - auto tracer_pkg = pmb->packages.Get("tracers"); - if (tracer_pkg->Param("enabled")) { - auto &mbd = pmb->meshblock_data.Get(); - Tracers::FillTracers(mbd.get()); - } } } // namespace turbulence diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 0cf5f204..5886bf65 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -20,6 +20,7 @@ #include "tracers.hpp" #include "../main.hpp" +#include "basic_types.hpp" #include "utils/error_checking.hpp" namespace Tracers { @@ -126,7 +127,7 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { * Registered Quantities (in addition to t, x, y, z): * rho, vel, B **/ -void FillTracers(MeshBlockData *mbd) { +TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { auto *pmb = mbd->GetParentPointer(); auto hydro_pkg = pmb->packages.Get("Hydro"); @@ -184,6 +185,7 @@ void FillTracers(MeshBlockData *mbd) { swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } }); + return TaskStatus::complete; } // FillTracers diff --git a/src/tracers/tracers.hpp b/src/tracers/tracers.hpp index 1a21b2a7..4008221e 100644 --- a/src/tracers/tracers.hpp +++ b/src/tracers/tracers.hpp @@ -39,7 +39,7 @@ std::shared_ptr Initialize(ParameterInput *pin); TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt); -void FillTracers(MeshBlockData *rc); +TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm); } // namespace Tracers From 0e86320b813dc2bcb8d21de99ac5520fa89223f0 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 6 Mar 2024 23:08:31 +0100 Subject: [PATCH 07/30] Added brief instructions for new particle field --- src/tracers/tracers.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 5886bf65..9e5fad63 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -22,6 +22,7 @@ #include "../main.hpp" #include "basic_types.hpp" #include "utils/error_checking.hpp" +#include namespace Tracers { using namespace parthenon::package::prelude; @@ -78,6 +79,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { tracer_pkg->AddSwarmValue("B_z", swarm_name, real_swarmvalue_metadata); } + // MARCUS AND EVAN LOOK + tracer_pkg->AddSwarmValue("lnrho_0", swarm_name, real_swarmvalue_metadata); + return tracer_pkg; } // Initialize @@ -155,6 +159,9 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { } auto &rho = swarm->Get("rho").Get(); auto &pressure = swarm->Get("pressure").Get(); + // MARCUS AND EVAN LOOK + const auto current_cycle = tm.ncycle; + auto &lnrho_0 = swarm->Get("lnrho_0").Get(); // Get hydro/mhd fluid vars const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); @@ -181,6 +188,11 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { B_z(n) = prim_pack(IB3, k, j, i); } + // MARCUS AND EVAN LOOK + if (current_cycle % 1 == 0) { + lnrho_0(n) = Kokkos::log(prim_pack(IDN, k, j, i)); + } + bool unsed_tmp = true; swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } From 4973d8df618b38651bfe53493975ec9682c81368 Mon Sep 17 00:00:00 2001 From: mbruggen Date: Wed, 6 Mar 2024 15:12:47 -0800 Subject: [PATCH 08/30] Update tracers.cpp filled in values for s_i --- src/tracers/tracers.cpp | 111 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 104 insertions(+), 7 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 9e5fad63..45bcc012 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -80,8 +80,31 @@ std::shared_ptr Initialize(ParameterInput *pin) { } // MARCUS AND EVAN LOOK - tracer_pkg->AddSwarmValue("lnrho_0", swarm_name, real_swarmvalue_metadata); - + tracer_pkg->AddSwarmValue("s_0", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_1", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_2", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_4", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_8", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_16", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_32", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_64", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_128", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_256", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_512", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("s_1024", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_0", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_1", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_2", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_4", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_8", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_16", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_32", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_64", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_128", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_256", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_512", swarm_name, real_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot_1024", swarm_name, real_swarmvalue_metadata); + return tracer_pkg; } // Initialize @@ -161,8 +184,32 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { auto &pressure = swarm->Get("pressure").Get(); // MARCUS AND EVAN LOOK const auto current_cycle = tm.ncycle; - auto &lnrho_0 = swarm->Get("lnrho_0").Get(); - + auto &s_0 = swarm->Get("s_0").Get(); + auto &s_1 = swarm->Get("s_1").Get(); + auto &s_2 = swarm->Get("s_2").Get(); + auto &s_4 = swarm->Get("s_4").Get(); + auto &s_8 = swarm->Get("s_8").Get(); + auto &s_16 = swarm->Get("s_16").Get(); + auto &s_32 = swarm->Get("s_32").Get(); + auto &s_64 = swarm->Get("s_64").Get(); + auto &s_128 = swarm->Get("s_128").Get(); + auto &s_256 = swarm->Get("s_256").Get(); + auto &s_512 = swarm->Get("s_512").Get(); + auto &s_1024 = swarm->Get("s_1024").Get(); + + auto &sdot_0 = swarm->Get("sdot_0").Get(); + auto &sdot_1 = swarm->Get("sdot_1").Get(); + auto &sdot_2 = swarm->Get("sdot_2").Get(); + auto &sdot_4 = swarm->Get("sdot_4").Get(); + auto &sdot_8 = swarm->Get("sdot_8").Get(); + auto &sdot_16 = swarm->Get("sdot_16").Get(); + auto &sdot_32 = swarm->Get("sdot_32").Get(); + auto &sdot_64 = swarm->Get("sdot_64").Get(); + auto &sdot_128 = swarm->Get("sdot_128").Get(); + auto &sdot_256 = swarm->Get("sdot_256").Get(); + auto &sdot_512 = swarm->Get("sdot_512").Get(); + auto &sdot_1024 = swarm->Get("sdot_1024").Get(); + // Get hydro/mhd fluid vars const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); @@ -188,11 +235,61 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { B_z(n) = prim_pack(IB3, k, j, i); } - // MARCUS AND EVAN LOOK + // MARCUS AND EVAN LOOK + // Q: DO WE HAVE TO INITIALISE S_0? + + if (current_cycle % 1024 == 0) { + s_1024(n) = s_512(n) + sdot_1024(n) = sdot_512(n) + + } + if (current_cycle % 512 == 0) { + s_512(n) = s_256(n) + sdot_512(n) = sdot_256(n) + } + if (current_cycle % 256 == 0) { + s_256(n) = s_128(n) + sdot_256(n) = sdot_128(n) + } + if (current_cycle % 128 == 0) { + s_128(n) = s_64(n) + sdot_128(n) = sdot_64(n) + + } + if (current_cycle % 64 == 0) { + s_64(n) = s_32(n) + sdot_64(n) = sdot_32(n) + } + if (current_cycle % 32 == 0) { + s_32(n) = s_16(n) + sdot_32(n) = sdot_16(n) + } + if (current_cycle % 16 == 0) { + s_16(n) = s_8(n) + sdot_16(n) = sdot_8(n) + } + if (current_cycle % 8 == 0) { + s_8(n) = s_4(n) + sdot_8(n) = sdot_4(n) + } + if (current_cycle % 4 == 0) { + s_4(n) = s_2(n) + sdot_4(n) = sdot_2(n) + } + if (current_cycle % 2 == 0) { + s_2(n) = s_1(n) + sdot_2(n) = sdot_1(n) + + // t2 = t1 here we build up our arrays of the previous times + } if (current_cycle % 1 == 0) { - lnrho_0(n) = Kokkos::log(prim_pack(IDN, k, j, i)); + s_1(n) = s_0(n) + sdot_1(n) = sdot_0(n) + // t1 = t0 here we build up our arrays of the previous times } - + s_0(n) = Kokkos::log(prim_pack(IDN, k, j, i)); + sdot_0(n) = (s_0(n)-s_1(n))/ DELTAT ; + bool unsed_tmp = true; swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } From 7c57381ca93469cffb7c59790a405ca3e423f971 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 7 Mar 2024 11:47:13 +0100 Subject: [PATCH 09/30] Simplify/dedup code --- src/tracers/tracers.cpp | 132 ++++++++-------------------------------- 1 file changed, 27 insertions(+), 105 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 45bcc012..258ec7cc 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -21,6 +21,7 @@ #include "tracers.hpp" #include "../main.hpp" #include "basic_types.hpp" +#include "interface/metadata.hpp" #include "utils/error_checking.hpp" #include @@ -49,7 +50,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Add swarm of tracers std::string swarm_name = "tracers"; - Metadata swarm_metadata({Metadata::Provides, Metadata::None}); + Metadata swarm_metadata({Metadata::Provides, Metadata::None, Metadata::Restart}); tracer_pkg->AddSwarm(swarm_name, swarm_metadata); Metadata real_swarmvalue_metadata({Metadata::Real}); tracer_pkg->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); @@ -80,31 +81,12 @@ std::shared_ptr Initialize(ParameterInput *pin) { } // MARCUS AND EVAN LOOK - tracer_pkg->AddSwarmValue("s_0", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_1", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_2", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_4", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_8", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_16", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_32", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_64", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_128", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_256", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_512", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("s_1024", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_0", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_1", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_2", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_4", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_8", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_16", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_32", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_64", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_128", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_256", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_512", swarm_name, real_swarmvalue_metadata); - tracer_pkg->AddSwarmValue("sdot_1024", swarm_name, real_swarmvalue_metadata); - + // Using a vector to reduce code duplication. + Metadata vreal_swarmvalue_metadata( + {Metadata::Real, Metadata::Vector, Metadata::Restart}, std::vector{12}); + tracer_pkg->AddSwarmValue("s", swarm_name, vreal_swarmvalue_metadata); + tracer_pkg->AddSwarmValue("sdot", swarm_name, vreal_swarmvalue_metadata); + return tracer_pkg; } // Initialize @@ -184,32 +166,10 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { auto &pressure = swarm->Get("pressure").Get(); // MARCUS AND EVAN LOOK const auto current_cycle = tm.ncycle; - auto &s_0 = swarm->Get("s_0").Get(); - auto &s_1 = swarm->Get("s_1").Get(); - auto &s_2 = swarm->Get("s_2").Get(); - auto &s_4 = swarm->Get("s_4").Get(); - auto &s_8 = swarm->Get("s_8").Get(); - auto &s_16 = swarm->Get("s_16").Get(); - auto &s_32 = swarm->Get("s_32").Get(); - auto &s_64 = swarm->Get("s_64").Get(); - auto &s_128 = swarm->Get("s_128").Get(); - auto &s_256 = swarm->Get("s_256").Get(); - auto &s_512 = swarm->Get("s_512").Get(); - auto &s_1024 = swarm->Get("s_1024").Get(); - - auto &sdot_0 = swarm->Get("sdot_0").Get(); - auto &sdot_1 = swarm->Get("sdot_1").Get(); - auto &sdot_2 = swarm->Get("sdot_2").Get(); - auto &sdot_4 = swarm->Get("sdot_4").Get(); - auto &sdot_8 = swarm->Get("sdot_8").Get(); - auto &sdot_16 = swarm->Get("sdot_16").Get(); - auto &sdot_32 = swarm->Get("sdot_32").Get(); - auto &sdot_64 = swarm->Get("sdot_64").Get(); - auto &sdot_128 = swarm->Get("sdot_128").Get(); - auto &sdot_256 = swarm->Get("sdot_256").Get(); - auto &sdot_512 = swarm->Get("sdot_512").Get(); - auto &sdot_1024 = swarm->Get("sdot_1024").Get(); - + const auto dt = tm.dt; + auto &s = swarm->Get("s").Get(); + auto &sdot = swarm->Get("sdot").Get(); + // Get hydro/mhd fluid vars const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); @@ -235,61 +195,23 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { B_z(n) = prim_pack(IB3, k, j, i); } - // MARCUS AND EVAN LOOK + // MARCUS AND EVAN LOOK // Q: DO WE HAVE TO INITIALISE S_0? - - if (current_cycle % 1024 == 0) { - s_1024(n) = s_512(n) - sdot_1024(n) = sdot_512(n) - - } - if (current_cycle % 512 == 0) { - s_512(n) = s_256(n) - sdot_512(n) = sdot_256(n) - } - if (current_cycle % 256 == 0) { - s_256(n) = s_128(n) - sdot_256(n) = sdot_128(n) - } - if (current_cycle % 128 == 0) { - s_128(n) = s_64(n) - sdot_128(n) = sdot_64(n) - - } - if (current_cycle % 64 == 0) { - s_64(n) = s_32(n) - sdot_64(n) = sdot_32(n) - } - if (current_cycle % 32 == 0) { - s_32(n) = s_16(n) - sdot_32(n) = sdot_16(n) - } - if (current_cycle % 16 == 0) { - s_16(n) = s_8(n) - sdot_16(n) = sdot_8(n) - } - if (current_cycle % 8 == 0) { - s_8(n) = s_4(n) - sdot_8(n) = sdot_4(n) - } - if (current_cycle % 4 == 0) { - s_4(n) = s_2(n) - sdot_4(n) = sdot_2(n) - } - if (current_cycle % 2 == 0) { - s_2(n) = s_1(n) - sdot_2(n) = sdot_1(n) - - // t2 = t1 here we build up our arrays of the previous times + // PG: It's default intiialized to zero, so probably not. + + int dncycle = 1024; + int s_idx = 11; + while (dncycle > 0) { + if (current_cycle % dncycle == 0) { + s(s_idx, n) = s(s_idx - 1, n); + sdot(s_idx, n) = sdot(s_idx - 1, n); + } + dncycle /= 2; + s_idx -= 1; } - if (current_cycle % 1 == 0) { - s_1(n) = s_0(n) - sdot_1(n) = sdot_0(n) - // t1 = t0 here we build up our arrays of the previous times - } - s_0(n) = Kokkos::log(prim_pack(IDN, k, j, i)); - sdot_0(n) = (s_0(n)-s_1(n))/ DELTAT ; - + s(0, n) = Kokkos::log(prim_pack(IDN, k, j, i)); + sdot(0, n) = (s(0, n) - s(1, n)) / dt; + bool unsed_tmp = true; swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } From 46b48e19e28aaa89f9bf4d978573785b8ff574c8 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 7 Mar 2024 15:08:19 +0100 Subject: [PATCH 10/30] Add lookback time --- src/tracers/tracers.cpp | 54 +++++++++++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 258ec7cc..ec8f2d54 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -18,12 +18,15 @@ // distribute copies to the public, perform publicly and display // publicly, and to permit others to do so. -#include "tracers.hpp" +#include +#include +#include + #include "../main.hpp" #include "basic_types.hpp" #include "interface/metadata.hpp" +#include "tracers.hpp" #include "utils/error_checking.hpp" -#include namespace Tracers { using namespace parthenon::package::prelude; @@ -50,10 +53,13 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Add swarm of tracers std::string swarm_name = "tracers"; + // TODO(pgrete) Check where metadata, e.g., for restart is required (i.e., at the swarm + // or variable level). Metadata swarm_metadata({Metadata::Provides, Metadata::None, Metadata::Restart}); tracer_pkg->AddSwarm(swarm_name, swarm_metadata); Metadata real_swarmvalue_metadata({Metadata::Real}); - tracer_pkg->AddSwarmValue("id", swarm_name, Metadata({Metadata::Integer})); + tracer_pkg->AddSwarmValue("id", swarm_name, + Metadata({Metadata::Integer, Metadata::Restart})); // TODO(pgrete) Add CheckDesired/required for vars // thermo variables @@ -80,12 +86,22 @@ std::shared_ptr Initialize(ParameterInput *pin) { tracer_pkg->AddSwarmValue("B_z", swarm_name, real_swarmvalue_metadata); } + // TODO(pgrete) this should eventually moved to a more pgen/user specific place // MARCUS AND EVAN LOOK + // Number of lookback times to be stored (in powers of 2, + // i.e., 12 allows to go from 0, 2^0 = 1, 2^1 = 2, 2^2 = 4, ..., 2^10 = 1024 cycles) + const int n_lookback = 12; // could even be made an input parameter if required/desired + // (though it should probably not be changeable for restarts) + tracer_pkg->AddParam("n_lookback", n_lookback); // Using a vector to reduce code duplication. Metadata vreal_swarmvalue_metadata( - {Metadata::Real, Metadata::Vector, Metadata::Restart}, std::vector{12}); + {Metadata::Real, Metadata::Vector, Metadata::Restart}, + std::vector{n_lookback}); tracer_pkg->AddSwarmValue("s", swarm_name, vreal_swarmvalue_metadata); tracer_pkg->AddSwarmValue("sdot", swarm_name, vreal_swarmvalue_metadata); + // Timestamps for the lookback entries + tracer_pkg->AddParam<>("t_lookback", std::vector(n_lookback), + Params::Mutability::Restart); return tracer_pkg; } // Initialize @@ -139,10 +155,10 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { auto *pmb = mbd->GetParentPointer(); - auto hydro_pkg = pmb->packages.Get("Hydro"); auto &sd = pmb->swarm_data.Get(); auto &swarm = sd->Get("tracers"); + auto hydro_pkg = pmb->packages.Get("Hydro"); const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) @@ -167,9 +183,32 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { // MARCUS AND EVAN LOOK const auto current_cycle = tm.ncycle; const auto dt = tm.dt; + auto &s = swarm->Get("s").Get(); auto &sdot = swarm->Get("sdot").Get(); + auto tracers_pkg = pmb->packages.Get("tracers"); + const auto n_lookback = tracers_pkg->Param("n_lookback"); + // Params (which is storing t_lookback) is shared across all blocks. Thus, we make sure + // to just call it once (for the first block with global id 0). + if (pmb->gid == 0) { + // Note, that this is a standard vector, so it cannot be used in the kernel (but also + // don't need to be used as can directly update it) + auto t_lookback = tracers_pkg->Param>("t_lookback"); + auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); + auto idx = n_lookback - 1; + while (dncycle > 0) { + if (current_cycle % dncycle == 0) { + t_lookback[idx] = t_lookback[idx - 1]; + } + dncycle /= 2; + idx -= 1; + } + t_lookback[0] = tm.time; + // Write data back to Params dict + tracers_pkg->UpdateParam("t_lookback", t_lookback); + } + // Get hydro/mhd fluid vars const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); @@ -198,9 +237,8 @@ TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { // MARCUS AND EVAN LOOK // Q: DO WE HAVE TO INITIALISE S_0? // PG: It's default intiialized to zero, so probably not. - - int dncycle = 1024; - int s_idx = 11; + auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); + auto s_idx = n_lookback - 1; while (dncycle > 0) { if (current_cycle % dncycle == 0) { s(s_idx, n) = s(s_idx - 1, n); From ddc5c789d6db3210ce612784c78164e333ec3b79 Mon Sep 17 00:00:00 2001 From: evan1022 Date: Fri, 8 Mar 2024 12:53:44 -0800 Subject: [PATCH 11/30] Update turbulence.cpp --- src/pgen/turbulence.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 23f101ee..f8138f39 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -536,5 +536,34 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto state_dist = few_modes_ft.GetDistState(); pin->SetString("problem/turbulence", "state_dist", state_dist); } +\\ +\\ COMPUTE THE AVERAGE OVER PRODUCTS OF S VALUES +\\ PHILIPP: PLEASE FIX THIS +\\ +auto correlation_s = tracers_pkg->Param>("correlation_s"); +auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); + +Kokkos::parallel_reduce( + "", + Kokkos::MDRangePolicy>(n), + KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { + auto tracers_pkg = pmb->packages.Get("tracers"); + const auto n_lookback = tracers_pkg->Param("n_lookback"); + auto idx = n_lookback - 1; + auto &s = swarm->Get("s").Get(); + auto &sdot = swarm->Get("sdot").Get(); + auto corr_s; + auto corr_sdot; + + for (i = 0; i < idx; i++) { + corr_s[i] =s(0, n) *s(i,n); \\PHILIPP: HOW DO YOU DO THIS INDEXING? + corr_sdot[i]=sdot(0, n)*sdot(i,n); + } + }, + Kokkos::Add(correlation_s,correlation_sdot)); + correlation_s=correlation_s/N; + correlation_sdot=correlation_sdot/N; + + } // namespace turbulence From 94c1eca60109a84ae90714d38612fdbe0e15ca4c Mon Sep 17 00:00:00 2001 From: evan1022 Date: Fri, 8 Mar 2024 12:56:13 -0800 Subject: [PATCH 12/30] Update turbulence.cpp --- src/pgen/turbulence.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index f8138f39..e4a36631 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -542,10 +542,10 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { \\ auto correlation_s = tracers_pkg->Param>("correlation_s"); auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); - Kokkos::parallel_reduce( "", Kokkos::MDRangePolicy>(n), + \\ PHILIPP Here is our attempt to create a kokkos lambda that returns the correlations of s and sdot for a single particle KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { auto tracers_pkg = pmb->packages.Get("tracers"); const auto n_lookback = tracers_pkg->Param("n_lookback"); @@ -560,6 +560,7 @@ Kokkos::parallel_reduce( corr_sdot[i]=sdot(0, n)*sdot(i,n); } }, + \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average Kokkos::Add(correlation_s,correlation_sdot)); correlation_s=correlation_s/N; correlation_sdot=correlation_sdot/N; From b81eb668eec9180f85f81d587f982c8716cf4208 Mon Sep 17 00:00:00 2001 From: mbruggen Date: Fri, 8 Mar 2024 14:35:21 -0800 Subject: [PATCH 13/30] Update turbulence.cpp small edits in the par_reduce in the computation of the correlations --- src/pgen/turbulence.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index e4a36631..eeba24b0 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -542,8 +542,9 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { \\ auto correlation_s = tracers_pkg->Param>("correlation_s"); auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); + Kokkos::parallel_reduce( - "", + "Correlation", Kokkos::MDRangePolicy>(n), \\ PHILIPP Here is our attempt to create a kokkos lambda that returns the correlations of s and sdot for a single particle KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { @@ -558,12 +559,11 @@ Kokkos::parallel_reduce( for (i = 0; i < idx; i++) { corr_s[i] =s(0, n) *s(i,n); \\PHILIPP: HOW DO YOU DO THIS INDEXING? corr_sdot[i]=sdot(0, n)*sdot(i,n); - } }, - \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average - Kokkos::Add(correlation_s,correlation_sdot)); - correlation_s=correlation_s/N; - correlation_sdot=correlation_sdot/N; + \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average + \\ N is the number of particles + Kokkos::Add(correlation_s,correlation_sdot); + correlation_s=correlation_s/N, correlation_sdot=correlation_sdot/N;) From 7be5d87a0878660718349397c2af5490238e3d6d Mon Sep 17 00:00:00 2001 From: mbruggen Date: Fri, 8 Mar 2024 14:41:25 -0800 Subject: [PATCH 14/30] Update tracers.cpp --- src/tracers/tracers.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index ec8f2d54..76dd848d 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -102,6 +102,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Timestamps for the lookback entries tracer_pkg->AddParam<>("t_lookback", std::vector(n_lookback), Params::Mutability::Restart); + tracer_pkg->AddParam<>("correlation_s", std::vector(n_lookback), + Params::Mutability::Restart); + tracer_pkg->AddParam<>("correlation_sdot", std::vector(n_lookback), + Params::Mutability::Restart); return tracer_pkg; } // Initialize From bc39143ce6787b79a29e0735569514793459fbd9 Mon Sep 17 00:00:00 2001 From: mbruggen Date: Fri, 8 Mar 2024 14:43:10 -0800 Subject: [PATCH 15/30] Update turbulence.cpp --- src/pgen/turbulence.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index eeba24b0..6b0ea193 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -563,8 +563,9 @@ Kokkos::parallel_reduce( \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average \\ N is the number of particles Kokkos::Add(correlation_s,correlation_sdot); - correlation_s=correlation_s/N, correlation_sdot=correlation_sdot/N;) - + correlation_s=correlation_s, correlation_sdot=correlation_sdot) + tracers_pkg->UpdateParam("correlation_s", correlation_s); + tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot); } // namespace turbulence From 582dc2aa7b8cc1177804115c6a4ea0b2ed3c0d7a Mon Sep 17 00:00:00 2001 From: mbruggen Date: Fri, 8 Mar 2024 14:46:19 -0800 Subject: [PATCH 16/30] Update turbulence.cpp --- src/pgen/turbulence.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 6b0ea193..3393279b 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -565,7 +565,7 @@ Kokkos::parallel_reduce( Kokkos::Add(correlation_s,correlation_sdot); correlation_s=correlation_s, correlation_sdot=correlation_sdot) - tracers_pkg->UpdateParam("correlation_s", correlation_s); - tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot); + tracers_pkg->UpdateParam("correlation_s", correlation_s/N); + tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot/N); } // namespace turbulence From 0a1a76ec076d4097212d356e5541c9a7c47f42f4 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 13:28:08 +0100 Subject: [PATCH 17/30] Use MeshData in FillTracers to allow reduction --- src/hydro/hydro_driver.cpp | 14 ++- src/tracers/tracers.cpp | 193 +++++++++++++++++++------------------ src/tracers/tracers.hpp | 2 +- 3 files changed, 110 insertions(+), 99 deletions(-) diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index b352c382..50b3385a 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -24,6 +24,7 @@ #include "glmmhd/glmmhd.hpp" #include "hydro.hpp" #include "hydro_driver.hpp" +#include "utils/error_checking.hpp" using namespace parthenon::driver::prelude; @@ -412,8 +413,17 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { auto receive = tl.AddTask(send, &SwarmContainer::Receive, sd.get(), BoundaryCommSubset::all); - - auto fill = tl.AddTask(receive, Tracers::FillTracers, mbd0.get(), tm); + } + // TODO(pgrete) Fix/cleanup once we got swarm packs. + // We need just a single region with a single task in order to be able to use plain + // reductions. + PARTHENON_REQUIRE_THROWS(num_partitions == 1, + "Only pack_size=-1 currently supported for tracers.") + TaskRegion &single_tasklist_per_pack_region_4 = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; i++) { + auto &tl = single_tasklist_per_pack_region_4[i]; + auto &mu0 = pmesh->mesh_data.GetOrAdd("base", i); + auto fill = tl.AddTask(none, Tracers::FillTracers, mu0.get(), tm); } } diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 76dd848d..56a265a4 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -156,108 +156,109 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { * Registered Quantities (in addition to t, x, y, z): * rho, vel, B **/ -TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm) { - - auto *pmb = mbd->GetParentPointer(); - auto &sd = pmb->swarm_data.Get(); - auto &swarm = sd->Get("tracers"); - - auto hydro_pkg = pmb->packages.Get("Hydro"); - const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; - - // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) - // pull swarm vars - auto &x = swarm->Get("x").Get(); - auto &y = swarm->Get("y").Get(); - auto &z = swarm->Get("z").Get(); - auto &vel_x = swarm->Get("vel_x").Get(); - auto &vel_y = swarm->Get("vel_y").Get(); - auto &vel_z = swarm->Get("vel_z").Get(); - // Assign some (definitely existing) default var - auto B_x = vel_x.Get(); - auto B_y = vel_x.Get(); - auto B_z = vel_x.Get(); - if (mhd) { - B_x = swarm->Get("B_x").Get(); - B_y = swarm->Get("B_y").Get(); - B_z = swarm->Get("B_z").Get(); - } - auto &rho = swarm->Get("rho").Get(); - auto &pressure = swarm->Get("pressure").Get(); - // MARCUS AND EVAN LOOK - const auto current_cycle = tm.ncycle; - const auto dt = tm.dt; - - auto &s = swarm->Get("s").Get(); - auto &sdot = swarm->Get("sdot").Get(); - - auto tracers_pkg = pmb->packages.Get("tracers"); - const auto n_lookback = tracers_pkg->Param("n_lookback"); - // Params (which is storing t_lookback) is shared across all blocks. Thus, we make sure - // to just call it once (for the first block with global id 0). - if (pmb->gid == 0) { - // Note, that this is a standard vector, so it cannot be used in the kernel (but also - // don't need to be used as can directly update it) - auto t_lookback = tracers_pkg->Param>("t_lookback"); - auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); - auto idx = n_lookback - 1; - while (dncycle > 0) { - if (current_cycle % dncycle == 0) { - t_lookback[idx] = t_lookback[idx - 1]; +TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { + // Get hydro/mhd fluid vars + const auto &prim_pack = md->PackVariables(std::vector{"prim"}); + for (int b = 0; b < md->NumBlocks(); b++) { + auto *pmb = md->GetBlockData(b)->GetBlockPointer(); + auto &sd = pmb->swarm_data.Get(); + auto &swarm = sd->Get("tracers"); + + auto hydro_pkg = pmb->packages.Get("Hydro"); + const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; + + // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) + // pull swarm vars + auto &x = swarm->Get("x").Get(); + auto &y = swarm->Get("y").Get(); + auto &z = swarm->Get("z").Get(); + auto &vel_x = swarm->Get("vel_x").Get(); + auto &vel_y = swarm->Get("vel_y").Get(); + auto &vel_z = swarm->Get("vel_z").Get(); + // Assign some (definitely existing) default var + auto B_x = vel_x.Get(); + auto B_y = vel_x.Get(); + auto B_z = vel_x.Get(); + if (mhd) { + B_x = swarm->Get("B_x").Get(); + B_y = swarm->Get("B_y").Get(); + B_z = swarm->Get("B_z").Get(); + } + auto &rho = swarm->Get("rho").Get(); + auto &pressure = swarm->Get("pressure").Get(); + // MARCUS AND EVAN LOOK + const auto current_cycle = tm.ncycle; + const auto dt = tm.dt; + + auto &s = swarm->Get("s").Get(); + auto &sdot = swarm->Get("sdot").Get(); + + auto tracers_pkg = pmb->packages.Get("tracers"); + const auto n_lookback = tracers_pkg->Param("n_lookback"); + // Params (which is storing t_lookback) is shared across all blocks. Thus, we make + // sure to just call it once (for the first block with global id 0). + if (pmb->gid == 0) { + // Note, that this is a standard vector, so it cannot be used in the kernel (but + // also don't need to be used as can directly update it) + auto t_lookback = tracers_pkg->Param>("t_lookback"); + auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); + auto idx = n_lookback - 1; + while (dncycle > 0) { + if (current_cycle % dncycle == 0) { + t_lookback[idx] = t_lookback[idx - 1]; + } + dncycle /= 2; + idx -= 1; } - dncycle /= 2; - idx -= 1; + t_lookback[0] = tm.time; + // Write data back to Params dict + tracers_pkg->UpdateParam("t_lookback", t_lookback); } - t_lookback[0] = tm.time; - // Write data back to Params dict - tracers_pkg->UpdateParam("t_lookback", t_lookback); - } - - // Get hydro/mhd fluid vars - const auto &prim_pack = mbd->PackVariables(std::vector{"prim"}); - auto swarm_d = swarm->GetDeviceContext(); - - // update loop. - const int max_active_index = swarm->GetMaxActiveIndex(); - pmb->par_for( - "Fill Tracers", 0, max_active_index, KOKKOS_LAMBDA(const int n) { - if (swarm_d.IsActive(n)) { - int k, j, i; - swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); - - // TODO(pgrete) Interpolate - rho(n) = prim_pack(IDN, k, j, i); - vel_x(n) = prim_pack(IV1, k, j, i); - vel_y(n) = prim_pack(IV2, k, j, i); - vel_z(n) = prim_pack(IV3, k, j, i); - pressure(n) = prim_pack(IPR, k, j, i); - if (mhd) { - B_x(n) = prim_pack(IB1, k, j, i); - B_y(n) = prim_pack(IB2, k, j, i); - B_z(n) = prim_pack(IB3, k, j, i); - } + auto swarm_d = swarm->GetDeviceContext(); + + // update loop. + const int max_active_index = swarm->GetMaxActiveIndex(); + pmb->par_for( + "Fill Tracers", 0, max_active_index, KOKKOS_LAMBDA(const int n) { + const auto prim = prim_pack(b); + if (swarm_d.IsActive(n)) { + int k, j, i; + swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); + + // TODO(pgrete) Interpolate + rho(n) = prim(IDN, k, j, i); + vel_x(n) = prim(IV1, k, j, i); + vel_y(n) = prim(IV2, k, j, i); + vel_z(n) = prim(IV3, k, j, i); + pressure(n) = prim(IPR, k, j, i); + if (mhd) { + B_x(n) = prim(IB1, k, j, i); + B_y(n) = prim(IB2, k, j, i); + B_z(n) = prim(IB3, k, j, i); + } - // MARCUS AND EVAN LOOK - // Q: DO WE HAVE TO INITIALISE S_0? - // PG: It's default intiialized to zero, so probably not. - auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); - auto s_idx = n_lookback - 1; - while (dncycle > 0) { - if (current_cycle % dncycle == 0) { - s(s_idx, n) = s(s_idx - 1, n); - sdot(s_idx, n) = sdot(s_idx - 1, n); + // MARCUS AND EVAN LOOK + // Q: DO WE HAVE TO INITIALISE S_0? + // PG: It's default intiialized to zero, so probably not. + auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); + auto s_idx = n_lookback - 1; + while (dncycle > 0) { + if (current_cycle % dncycle == 0) { + s(s_idx, n) = s(s_idx - 1, n); + sdot(s_idx, n) = sdot(s_idx - 1, n); + } + dncycle /= 2; + s_idx -= 1; } - dncycle /= 2; - s_idx -= 1; - } - s(0, n) = Kokkos::log(prim_pack(IDN, k, j, i)); - sdot(0, n) = (s(0, n) - s(1, n)) / dt; + s(0, n) = Kokkos::log(prim(IDN, k, j, i)); + sdot(0, n) = (s(0, n) - s(1, n)) / dt; - bool unsed_tmp = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); - } - }); + bool unsed_tmp = true; + swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); + } + }); + } return TaskStatus::complete; } // FillTracers diff --git a/src/tracers/tracers.hpp b/src/tracers/tracers.hpp index 4008221e..fe9c68d7 100644 --- a/src/tracers/tracers.hpp +++ b/src/tracers/tracers.hpp @@ -39,7 +39,7 @@ std::shared_ptr Initialize(ParameterInput *pin); TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt); -TaskStatus FillTracers(MeshBlockData *mbd, parthenon::SimTime &tm); +TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm); } // namespace Tracers From 62b759aaab349b4495822df4eae443fa3c415ff3 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 14:43:54 +0100 Subject: [PATCH 18/30] Update FillTracers internals to loop over blocks logic --- src/pgen/turbulence.cpp | 32 ------------- src/tracers/tracers.cpp | 100 +++++++++++++++++++++++++--------------- 2 files changed, 62 insertions(+), 70 deletions(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 3393279b..677186c2 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -536,36 +536,4 @@ void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) { auto state_dist = few_modes_ft.GetDistState(); pin->SetString("problem/turbulence", "state_dist", state_dist); } -\\ -\\ COMPUTE THE AVERAGE OVER PRODUCTS OF S VALUES -\\ PHILIPP: PLEASE FIX THIS -\\ -auto correlation_s = tracers_pkg->Param>("correlation_s"); -auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); - -Kokkos::parallel_reduce( - "Correlation", - Kokkos::MDRangePolicy>(n), - \\ PHILIPP Here is our attempt to create a kokkos lambda that returns the correlations of s and sdot for a single particle - KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { - auto tracers_pkg = pmb->packages.Get("tracers"); - const auto n_lookback = tracers_pkg->Param("n_lookback"); - auto idx = n_lookback - 1; - auto &s = swarm->Get("s").Get(); - auto &sdot = swarm->Get("sdot").Get(); - auto corr_s; - auto corr_sdot; - - for (i = 0; i < idx; i++) { - corr_s[i] =s(0, n) *s(i,n); \\PHILIPP: HOW DO YOU DO THIS INDEXING? - corr_sdot[i]=sdot(0, n)*sdot(i,n); - }, - \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average - \\ N is the number of particles - Kokkos::Add(correlation_s,correlation_sdot); - correlation_s=correlation_s, correlation_sdot=correlation_sdot) - - tracers_pkg->UpdateParam("correlation_s", correlation_s/N); - tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot/N); - } // namespace turbulence diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 56a265a4..26a46e2f 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -102,10 +102,6 @@ std::shared_ptr Initialize(ParameterInput *pin) { // Timestamps for the lookback entries tracer_pkg->AddParam<>("t_lookback", std::vector(n_lookback), Params::Mutability::Restart); - tracer_pkg->AddParam<>("correlation_s", std::vector(n_lookback), - Params::Mutability::Restart); - tracer_pkg->AddParam<>("correlation_sdot", std::vector(n_lookback), - Params::Mutability::Restart); return tracer_pkg; } // Initialize @@ -143,6 +139,9 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { y(n) += prim_pack(IV2, k, j, i) * dt; z(n) += prim_pack(IV3, k, j, i) * dt; + // The following call is required as it updates the internal block id following + // the advection. The internal id is used in the subsequent task to communicate + // particles. bool unused_temp = true; swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unused_temp); } @@ -157,16 +156,38 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { * rho, vel, B **/ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { - // Get hydro/mhd fluid vars + const auto current_cycle = tm.ncycle; + const auto dt = tm.dt; + + auto hydro_pkg = md->GetParentPointer()->packages.Get("Hydro"); + const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; + + auto tracers_pkg = md->GetParentPointer()->packages.Get("tracers"); + const auto n_lookback = tracers_pkg->Param("n_lookback"); + // Params (which is storing t_lookback) is shared across all blocks so we update it + // outside the block loop. Note, that this is a standard vector, so it cannot be used in + // the kernel (but also don't need to be used as can directly update it) + auto t_lookback = tracers_pkg->Param>("t_lookback"); + auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); + auto idx = n_lookback - 1; + while (dncycle > 0) { + if (current_cycle % dncycle == 0) { + t_lookback[idx] = t_lookback[idx - 1]; + } + dncycle /= 2; + idx -= 1; + } + t_lookback[0] = tm.time; + // Write data back to Params dict + tracers_pkg->UpdateParam("t_lookback", t_lookback); + + // Get hydro/mhd fluid vars over all blocks const auto &prim_pack = md->PackVariables(std::vector{"prim"}); for (int b = 0; b < md->NumBlocks(); b++) { auto *pmb = md->GetBlockData(b)->GetBlockPointer(); auto &sd = pmb->swarm_data.Get(); auto &swarm = sd->Get("tracers"); - auto hydro_pkg = pmb->packages.Get("Hydro"); - const auto mhd = hydro_pkg->Param("fluid") == Fluid::glmmhd; - // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) // pull swarm vars auto &x = swarm->Get("x").Get(); @@ -187,34 +208,9 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { auto &rho = swarm->Get("rho").Get(); auto &pressure = swarm->Get("pressure").Get(); // MARCUS AND EVAN LOOK - const auto current_cycle = tm.ncycle; - const auto dt = tm.dt; - auto &s = swarm->Get("s").Get(); auto &sdot = swarm->Get("sdot").Get(); - auto tracers_pkg = pmb->packages.Get("tracers"); - const auto n_lookback = tracers_pkg->Param("n_lookback"); - // Params (which is storing t_lookback) is shared across all blocks. Thus, we make - // sure to just call it once (for the first block with global id 0). - if (pmb->gid == 0) { - // Note, that this is a standard vector, so it cannot be used in the kernel (but - // also don't need to be used as can directly update it) - auto t_lookback = tracers_pkg->Param>("t_lookback"); - auto dncycle = static_cast(Kokkos::pow(2, n_lookback - 2)); - auto idx = n_lookback - 1; - while (dncycle > 0) { - if (current_cycle % dncycle == 0) { - t_lookback[idx] = t_lookback[idx - 1]; - } - dncycle /= 2; - idx -= 1; - } - t_lookback[0] = tm.time; - // Write data back to Params dict - tracers_pkg->UpdateParam("t_lookback", t_lookback); - } - auto swarm_d = swarm->GetDeviceContext(); // update loop. @@ -253,14 +249,42 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { } s(0, n) = Kokkos::log(prim(IDN, k, j, i)); sdot(0, n) = (s(0, n) - s(1, n)) / dt; - - bool unsed_tmp = true; - swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), unsed_tmp); } }); +#if 0 +\\ +\\ COMPUTE THE AVERAGE OVER PRODUCTS OF S VALUES +\\ PHILIPP: PLEASE FIX THIS +\\ +auto correlation_s = tracers_pkg->Param>("correlation_s"); +auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); + +Kokkos::parallel_reduce( + "Correlation", + Kokkos::MDRangePolicy>(n), + \\ PHILIPP Here is our attempt to create a kokkos lambda that returns the correlations of s and sdot for a single particle + KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { + auto tracers_pkg = pmb->packages.Get("tracers"); + const auto n_lookback = tracers_pkg->Param("n_lookback"); + auto idx = n_lookback - 1; + auto &s = swarm->Get("s").Get(); + auto &sdot = swarm->Get("sdot").Get(); + auto corr_s; + auto corr_sdot; + + for (i = 0; i < idx; i++) { + corr_s[i] =s(0, n) *s(i,n); \\PHILIPP: HOW DO YOU DO THIS INDEXING? + corr_sdot[i]=sdot(0, n)*sdot(i,n); + }, + \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average + \\ N is the number of particles + Kokkos::Add(correlation_s,correlation_sdot); + correlation_s=correlation_s, correlation_sdot=correlation_sdot) + + tracers_pkg->UpdateParam("correlation_s", correlation_s/N); + tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot/N); +#endif } return TaskStatus::complete; - } // FillTracers - } // namespace Tracers From 78ac98a1753cefcb695884229f82b0b1f5601d46 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 14:55:17 +0100 Subject: [PATCH 19/30] Calc correlations --- src/tracers/tracers.cpp | 49 +++++++++++++---------------------------- 1 file changed, 15 insertions(+), 34 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 26a46e2f..8982736a 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -181,6 +181,12 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // Write data back to Params dict tracers_pkg->UpdateParam("t_lookback", t_lookback); + // TODO(pgrete) Benchmark atomic and potentially update to proper reduction instead of + // atomics. + // Used for the parallel reduction. Could be reused but this way it's initalized to 0. + parthenon::ParArray2D corr("tracer correlations", 2, n_lookback); + int64_t num_particles_total = 0; + // Get hydro/mhd fluid vars over all blocks const auto &prim_pack = md->PackVariables(std::vector{"prim"}); for (int b = 0; b < md->NumBlocks(); b++) { @@ -207,7 +213,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { } auto &rho = swarm->Get("rho").Get(); auto &pressure = swarm->Get("pressure").Get(); - // MARCUS AND EVAN LOOK + auto &s = swarm->Get("s").Get(); auto &sdot = swarm->Get("sdot").Get(); @@ -249,41 +255,16 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { } s(0, n) = Kokkos::log(prim(IDN, k, j, i)); sdot(0, n) = (s(0, n) - s(1, n)) / dt; + + // Now that all s and sdot entries are updated, we calculate the (mean) + // correlations + for (s_idx = 0; s_idx < n_lookback; s_idx++) { + Kokkos::atomic_add(&corr(0, s_idx), s(0, n) * s(s_idx, n)); + Kokkos::atomic_add(&corr(1, s_idx), sdot(0, n) * s(s_idx, n)); + } } }); -#if 0 -\\ -\\ COMPUTE THE AVERAGE OVER PRODUCTS OF S VALUES -\\ PHILIPP: PLEASE FIX THIS -\\ -auto correlation_s = tracers_pkg->Param>("correlation_s"); -auto correlation_sdot = tracers_pkg->Param>("correlation_sdot"); - -Kokkos::parallel_reduce( - "Correlation", - Kokkos::MDRangePolicy>(n), - \\ PHILIPP Here is our attempt to create a kokkos lambda that returns the correlations of s and sdot for a single particle - KOKKOS_LAMBDA(const int n, Real &corr_s, Real &corr_sdot) { - auto tracers_pkg = pmb->packages.Get("tracers"); - const auto n_lookback = tracers_pkg->Param("n_lookback"); - auto idx = n_lookback - 1; - auto &s = swarm->Get("s").Get(); - auto &sdot = swarm->Get("sdot").Get(); - auto corr_s; - auto corr_sdot; - - for (i = 0; i < idx; i++) { - corr_s[i] =s(0, n) *s(i,n); \\PHILIPP: HOW DO YOU DO THIS INDEXING? - corr_sdot[i]=sdot(0, n)*sdot(i,n); - }, - \\ PHILLIP Here we are trying to do a reduction over all the particles that gives the average - \\ N is the number of particles - Kokkos::Add(correlation_s,correlation_sdot); - correlation_s=correlation_s, correlation_sdot=correlation_sdot) - - tracers_pkg->UpdateParam("correlation_s", correlation_s/N); - tracers_pkg->UpdateParam("correlation_sdot", correlation_sdot/N); -#endif + num_particles_total += swarm->GetNumActive(); } return TaskStatus::complete; } // FillTracers From b1ee6c029f5634abdf0f96fe9b9f27422ad07c1e Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 15:51:49 +0100 Subject: [PATCH 20/30] Dump corr to file --- src/tracers/tracers.cpp | 52 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 8982736a..65e53225 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -19,12 +19,15 @@ // publicly, and to permit others to do so. #include +#include #include #include #include "../main.hpp" #include "basic_types.hpp" #include "interface/metadata.hpp" +#include "kokkos_abstraction.hpp" +#include "parthenon_array_generic.hpp" #include "tracers.hpp" #include "utils/error_checking.hpp" @@ -265,7 +268,56 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { } }); num_particles_total += swarm->GetNumActive(); + } // loop over all blocks on this rank (this MeshData container) + + // Results still live in device memory. Copy to host for global reduction and output. + auto corr_h = Kokkos::create_mirror_view_and_copy(parthenon::HostMemSpace(), corr); +#ifdef MPI_PARALLEL + if (parthenon::Globals::my_rank == 0) { + PARTHENON_MPI_CHECK(MPI_Reduce(MPI_IN_PLACE, corr_h.data(), corr_h.GetSize(), + MPI_PARTHENON_REAL, MPI_SUM, 0, MPI_COMM_WORLD)); + PARTHENON_MPI_CHECK(MPI_Reduce(MPI_IN_PLACE, &num_particles_total, 1, MPI_INT64_T, + MPI_SUM, 0, MPI_COMM_WORLD)); + } else { + PARTHENON_MPI_CHECK(MPI_Reduce(corr_h.data(), corr_h.data(), corr_h.GetSize(), + MPI_PARTHENON_REAL, MPI_SUM, 0, MPI_COMM_WORLD)); + PARTHENON_MPI_CHECK(MPI_Reduce(&num_particles_total, &num_particles_total, 1, + MPI_INT64_T, MPI_SUM, 0, MPI_COMM_WORLD)); } +#endif + if (parthenon::Globals::my_rank == 0) { + // Turn sum into mean + for (int i = 0; i < n_lookback; i++) { + corr_h(0, i) /= static_cast(num_particles_total); + corr_h(1, i) /= static_cast(num_particles_total); + } + + // and write data + std::ofstream outfile; + const std::string fname("correlations.csv"); + // On startup, write header + if (current_cycle == 0) { + outfile.open(fname, std::ofstream::out); + outfile << "# Hello world\n"; + } else { + outfile.open(fname, std::ofstream::out | std::ofstream::app); + } + + outfile << tm.ncycle << "," << tm.time; + + for (int j = 0; j < 2; j++) { + for (int i = 0; i < n_lookback; i++) { + outfile << "," << corr_h(j, i); + } + } + for (int i = 0; i < n_lookback; i++) { + outfile << "," << t_lookback[i]; + } + outfile << std::endl; + + outfile.close(); + } + return TaskStatus::complete; } // FillTracers } // namespace Tracers From 4eeffab1e93efb796550d64f8e5311a9eb47dac4 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 20:24:45 +0100 Subject: [PATCH 21/30] Import interpolation from Phoebus --- src/CMakeLists.txt | 2 + src/tracers/tracers.cpp | 1 + src/utils/interpolation.hpp | 238 ++++++++++++++++++++++++++++++++++++ 3 files changed, 241 insertions(+) create mode 100644 src/utils/interpolation.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 40c5eb40..de101a6f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,6 +21,8 @@ add_executable( tracers/tracers.cpp tracers/tracers.hpp utils/few_modes_ft.cpp + utils/few_modes_ft.hpp + utils/interpolation.hpp ) add_subdirectory(pgen) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 65e53225..d4e00484 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -30,6 +30,7 @@ #include "parthenon_array_generic.hpp" #include "tracers.hpp" #include "utils/error_checking.hpp" +#include "../utils/interpolation.hpp" namespace Tracers { using namespace parthenon::package::prelude; diff --git a/src/utils/interpolation.hpp b/src/utils/interpolation.hpp new file mode 100644 index 00000000..8cf17927 --- /dev/null +++ b/src/utils/interpolation.hpp @@ -0,0 +1,238 @@ + +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Interpolation copied/refactored from +// https://github.com/lanl/phoebus and https://github.com/lanl/spiner +//======================================================================================== +// © 2022. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract +// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +// is operated by Triad National Security, LLC for the U.S. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#ifndef UTILS_INTERPOLATION_HPP_ +#define UTILS_INTERPOLATION_HPP_ + +// Spiner includes +#include + +// Parthenon includes +#include +#include +#include + +// Phoebus includes +#include "phoebus_utils/grid_utils.hpp" +#include "phoebus_utils/robust.hpp" + +namespace interpolation { + +using namespace parthenon::package::prelude; +using parthenon::Coordinates_t; +using weights_t = Spiner::weights_t; + +/// Base class for providing interpolation methods on uniformly spaced data. +/// Constructor is provided with spacing, number of support points, and desired +/// shift. GetIndicesAndWeights then updates arrays of indices and weights for +/// calculating the interpolated data. These arrays are of size StencilSize(). +/// Data is forced to zero outside the boundaries. +class Interpolation { + public: + KOKKOS_FUNCTION + Interpolation(const int nSupport, const Real dx, const Real shift) + : nSupport_(nSupport), dx_(dx), shift_(shift), ishift_(std::round(shift)) {} + + KOKKOS_INLINE_FUNCTION + virtual void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const {} + KOKKOS_INLINE_FUNCTION + virtual int StencilSize() const { return 0; } + + static constexpr int maxStencilSize = 2; + + protected: + const int nSupport_; + const Real dx_; + Real shift_; + int ishift_; +}; + +class PiecewiseConstant : public Interpolation { + public: + KOKKOS_FUNCTION + PiecewiseConstant(const int nSupport, const Real dx, const Real shift) + : Interpolation(nSupport, dx, shift) {} + + KOKKOS_INLINE_FUNCTION + void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { + idx[0] = i + ishift_; + wgt[0] = 1.; + if (idx[0] < 0 || idx[0] >= nSupport_) { + idx[0] = 0; + wgt[0] = 0.; + } + } + + KOKKOS_INLINE_FUNCTION + int StencilSize() const override { return 1; } +}; + +class Linear : public Interpolation { + public: + KOKKOS_FUNCTION + Linear(const int nSupport, const Real dx, const Real shift) + : Interpolation(nSupport, dx, shift) { + PARTHENON_FAIL("Not written yet!"); + } + + KOKKOS_INLINE_FUNCTION + void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { + idx[0] = std::floor(i + shift_); + idx[1] = idx[0] + 1; + + wgt[0] = wgt[1] = 1. - wgt[0]; + + for (int nsup = 0; nsup < 2; nsup++) { + if (idx[nsup] < 0 || idx[nsup] >= nSupport_) { + idx[nsup] = 0; + wgt[nsup] = 0.; + } + } + } + + KOKKOS_INLINE_FUNCTION + int StencilSize() const override { return 2; } +}; + +// TODO(JMM): Is this interpolation::Do syntax reasonable? An +// alternative path would be a class called "LCInterp with all +// static functions. Then it could have an `operator()` which would +// be maybe nicer? +// TODO(JMM): Merge this w/ what Ben has done. +namespace Cent { +namespace Linear { + +/* + * Get interpolation weights for linear interpolation + * PARAM[IN] - x - location to interpolate to + * PARAM[IN] - nx - number of points along this direction. Used for sanity checks. + * PARAM[IN] - coords - parthenon coords object + * PARAM[OUT] - ix - index of points to interpolate + * PARAM[OUT] - w - weights + */ +template +KOKKOS_INLINE_FUNCTION void GetWeights(const Real x, const int nx, + const Coordinates_t &coords, int &ix, + weights_t &w) { + const Real min = Coordinates::GetXv(0, coords); + const Real dx = coords.CellWidthFA(DIR); + ix = std::min(std::max(0, static_cast(robust::ratio(x - min, dx))), nx - 2); + const Real floor = min + ix * dx; + w[1] = robust::ratio(x - floor, dx); + w[0] = 1. - w[1]; +} + +/* + * Trilinear interpolation on a variable or meshblock pack + * PARAM[IN] - b - Meshblock index + * PARAM[IN] - X1, X2, X3 - Coordinate locations + * PARAM[IN] - p - Variable or MeshBlockPack + * PARAM[IN] - v - variable index + */ +template +KOKKOS_INLINE_FUNCTION Real Do3D(int b, const Real X1, const Real X2, const Real X3, + const Pack &p, int v) { + const auto &coords = p.GetCoords(b); + int ix[3]; + weights_t w[3]; + GetWeights(X1, p.GetDim(1), coords, ix[0], w[0]); + GetWeights(X2, p.GetDim(2), coords, ix[1], w[1]); + GetWeights(X3, p.GetDim(3), coords, ix[2], w[2]); + return (w[2][0] * (w[1][0] * (w[0][0] * p(b, v, ix[2], ix[1], ix[0]) + + w[0][1] * p(b, v, ix[2], ix[1], ix[0] + 1)) + + w[1][1] * (w[0][0] * p(b, v, ix[2], ix[1] + 1, ix[0]) + + w[0][1] * p(b, v, ix[2], ix[1] + 1, ix[0] + 1))) + + w[2][1] * (w[1][0] * (w[0][0] * p(b, v, ix[2] + 1, ix[1], ix[0]) + + w[0][1] * p(b, v, ix[2] + 1, ix[1], ix[0] + 1)) + + w[1][1] * (w[0][0] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0]) + + w[0][1] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0] + 1)))); +} + +/* + * Bilinear interpolation on a variable or meshblock pack + * PARAM[IN] - b - Meshblock index + * PARAM[IN] - X1, X2 - Coordinate locations + * PARAM[IN] - p - Variable or MeshBlockPack + * PARAM[IN] - v - variable index + */ +template +KOKKOS_INLINE_FUNCTION Real Do2D(int b, const Real X1, const Real X2, const Pack &p, + int v) { + const auto &coords = p.GetCoords(b); + int ix1, ix2; + weights_t w1, w2; + GetWeights(X1, p.GetDim(1), coords, ix1, w1); + GetWeights(X2, p.GetDim(2), coords, ix2, w2); + return (w2[0] * (w1[0] * p(b, v, 0, ix2, ix1) + w1[1] * p(b, v, 0, ix2, ix1 + 1)) + + w2[1] * + (w1[0] * p(b, v, 0, ix2 + 1, ix1) + w1[1] * p(b, v, 0, ix2 + 1, ix1 + 1))); +} + +/* + * Linear interpolation on a variable or meshblock pack + * PARAM[IN] - b - Meshblock index + * PARAM[IN] - X1 - Coordinate location + * PARAM[IN] - p - Variable or MeshBlockPack + * PARAM[IN] - v - variable index + */ +template +KOKKOS_INLINE_FUNCTION Real Do1D(int b, const Real X1, const Pack &p, int v) { + const auto &coords = p.GetCoords(b); + int ix; + weights_t w; + GetWeights(X1, p.GetDim(1), coords, ix, w); + return w[0] * p(b, v, 0, 0, ix) + w[1] * p(b, v, 0, 0, ix + 1); +} + +/* + * Trilinear or bilinear interpolation on a variable or meshblock pack + * PARAM[IN] - axisymmetric + * PARAM[IN] - b - Meshblock index + * PARAM[IN] - X1, X2, X3 - Coordinate locations + * PARAM[IN] - p - Variable or MeshBlockPack + * PARAM[IN] - v - variable index + */ +// JMM: I know this won't vectorize because of the switch, but it +// probably won't anyway, since we're doing trilinear +// interpolation, which will kill memory locality. Doing it this +// way means we can do trilinear vs bilinear which I think is a +// sufficient win at minimum code bloat. +template +KOKKOS_INLINE_FUNCTION Real Do(int b, const Real X1, const Real X2, const Real X3, + const Pack &p, int v) { + if (p.GetDim(3) > 1) { + return Do3D(b, X1, X2, X3, p, v); + } else if (p.GetDim(2) > 1) { + return Do2D(b, X1, X2, p, v); + } else { // 1D + return Do1D(b, X1, p, v); + } +} + +} // namespace Linear +} // namespace Cent +} // namespace interpolation + +// Convenience Namespace Alias +namespace LCInterp = interpolation::Cent::Linear; + +#endif // PHOEBUS_UTILS_PHOEBUS_INTERPOLATION_HPP_ From 59dbdfae51598acf9b72c5e5d7fb0e33ea5ea986 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 20:49:25 +0100 Subject: [PATCH 22/30] Adjust interface --- src/CMakeLists.txt | 1 + src/utils/interpolation.hpp | 27 +++++++++++++++++---------- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index de101a6f..05efed4d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,6 +23,7 @@ add_executable( utils/few_modes_ft.cpp utils/few_modes_ft.hpp utils/interpolation.hpp + utils/robust.hpp ) add_subdirectory(pgen) diff --git a/src/utils/interpolation.hpp b/src/utils/interpolation.hpp index 8cf17927..6185ffc9 100644 --- a/src/utils/interpolation.hpp +++ b/src/utils/interpolation.hpp @@ -1,4 +1,3 @@ - //======================================================================================== // AthenaPK - a performance portable block structured AMR astrophysical MHD code. // Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. @@ -23,23 +22,28 @@ #ifndef UTILS_INTERPOLATION_HPP_ #define UTILS_INTERPOLATION_HPP_ -// Spiner includes -#include - // Parthenon includes +#include "utils/error_checking.hpp" #include #include #include -// Phoebus includes -#include "phoebus_utils/grid_utils.hpp" -#include "phoebus_utils/robust.hpp" +#include "robust.hpp" namespace interpolation { using namespace parthenon::package::prelude; using parthenon::Coordinates_t; -using weights_t = Spiner::weights_t; + +// From https://github.com/lanl/spiner/blob/main/spiner/regular_grid_1d.hpp +// a poor-man's std::pair +struct weights_t { + Real first, second; + KOKKOS_INLINE_FUNCTION Real &operator[](const int i) { + assert(0 <= i && i <= 1); + return i == 0 ? first : second; + } +}; /// Base class for providing interpolation methods on uniformly spaced data. /// Constructor is provided with spacing, number of support points, and desired @@ -133,7 +137,10 @@ template KOKKOS_INLINE_FUNCTION void GetWeights(const Real x, const int nx, const Coordinates_t &coords, int &ix, weights_t &w) { - const Real min = Coordinates::GetXv(0, coords); + PARTHENON_DEBUG_REQUIRE( + typeid(Coordinates_t) == typeid(UniformCartesian), + "Interpolation routines currently only work for UniformCartesian"); + const Real min = coords.Xc(0); // assume uniform Cartesian const Real dx = coords.CellWidthFA(DIR); ix = std::min(std::max(0, static_cast(robust::ratio(x - min, dx))), nx - 2); const Real floor = min + ix * dx; @@ -235,4 +242,4 @@ KOKKOS_INLINE_FUNCTION Real Do(int b, const Real X1, const Real X2, const Real X // Convenience Namespace Alias namespace LCInterp = interpolation::Cent::Linear; -#endif // PHOEBUS_UTILS_PHOEBUS_INTERPOLATION_HPP_ +#endif // UTILS_INTERPOLATION_HPP_ From 71c2cd4d45a9a09f50da66c6b538ea74a8ea96e8 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 20:51:49 +0100 Subject: [PATCH 23/30] User interpolation for FillTracers --- src/tracers/tracers.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index d4e00484..b9f518de 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -24,13 +24,13 @@ #include #include "../main.hpp" +#include "../utils/interpolation.hpp" #include "basic_types.hpp" #include "interface/metadata.hpp" #include "kokkos_abstraction.hpp" #include "parthenon_array_generic.hpp" #include "tracers.hpp" #include "utils/error_checking.hpp" -#include "../utils/interpolation.hpp" namespace Tracers { using namespace parthenon::package::prelude; @@ -227,21 +227,20 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { const int max_active_index = swarm->GetMaxActiveIndex(); pmb->par_for( "Fill Tracers", 0, max_active_index, KOKKOS_LAMBDA(const int n) { - const auto prim = prim_pack(b); if (swarm_d.IsActive(n)) { int k, j, i; swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); // TODO(pgrete) Interpolate - rho(n) = prim(IDN, k, j, i); - vel_x(n) = prim(IV1, k, j, i); - vel_y(n) = prim(IV2, k, j, i); - vel_z(n) = prim(IV3, k, j, i); - pressure(n) = prim(IPR, k, j, i); + rho(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IDN); + vel_x(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IV1); + vel_y(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IV2); + vel_z(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IV3); + pressure(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IPR); if (mhd) { - B_x(n) = prim(IB1, k, j, i); - B_y(n) = prim(IB2, k, j, i); - B_z(n) = prim(IB3, k, j, i); + B_x(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IB1); + B_y(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IB2); + B_z(n) = LCInterp::Do(b, x(n), y(n), z(n), prim_pack, IB3); } // MARCUS AND EVAN LOOK @@ -257,7 +256,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { dncycle /= 2; s_idx -= 1; } - s(0, n) = Kokkos::log(prim(IDN, k, j, i)); + s(0, n) = Kokkos::log(rho(n)); sdot(0, n) = (s(0, n) - s(1, n)) / dt; // Now that all s and sdot entries are updated, we calculate the (mean) From 6bb07d7b3932c2905ed4482bb76b8c9df23154e0 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 11 Mar 2024 21:07:52 +0100 Subject: [PATCH 24/30] Use 2nd order time integrator for particles --- src/tracers/tracers.cpp | 33 ++++++++++++++------- src/utils/robust.hpp | 66 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 11 deletions(-) create mode 100644 src/utils/robust.hpp diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index b9f518de..cca35afd 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -118,6 +118,9 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { auto &x = swarm->Get("x").Get(); auto &y = swarm->Get("y").Get(); auto &z = swarm->Get("z").Get(); + auto &vel_x = swarm->Get("vel_x").Get(); + auto &vel_y = swarm->Get("vel_y").Get(); + auto &vel_z = swarm->Get("vel_z").Get(); auto swarm_d = swarm->GetDeviceContext(); @@ -131,17 +134,25 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { int k, j, i; swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k); - // Predictor/corrector will first make sense with non constant interpolation - // TODO(pgrete) add non-cnonst interpolation - // predictor - // const Real kx = x(n) + 0.5 * dt * rhs1; - // const Real ky = y(n) + 0.5 * dt * rhs2; - // const Real kz = z(n) + 0.5 * dt * rhs3; - - // corrector - x(n) += prim_pack(IV1, k, j, i) * dt; - y(n) += prim_pack(IV2, k, j, i) * dt; - z(n) += prim_pack(IV3, k, j, i) * dt; + // RK2/Heun's method (as the default in Flash) + // https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node130.html#SECTION06813000000000000000 + // Intermediate position and velocities + // x^{*,n+1} = x^n + dt * v^n + const auto x_star = x(n) + dt * vel_x(n); + const auto y_star = y(n) + dt * vel_y(n); + const auto z_star = z(n) + dt * vel_z(n); + + // v^{*,n+1} = v(x^{*,n+1}, t^{n+1}) + // First parameter b=0 assume to operate on a pack of a single block and needs + // to be updated if this becomes a MeshData function + const auto vel_x_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV1); + const auto vel_y_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV2); + const auto vel_z_star = LCInterp::Do(0, x_star, y_star, z_star, prim_pack, IV3); + + // Full update using mean velocity + x(n) += dt * 0.5 * (vel_x(n) + vel_x_star); + y(n) += dt * 0.5 * (vel_y(n) + vel_y_star); + z(n) += dt * 0.5 * (vel_z(n) + vel_z_star); // The following call is required as it updates the internal block id following // the advection. The internal id is used in the subsequent task to communicate diff --git a/src/utils/robust.hpp b/src/utils/robust.hpp new file mode 100644 index 00000000..7488abea --- /dev/null +++ b/src/utils/robust.hpp @@ -0,0 +1,66 @@ + +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Copied from https://github.com/lanl/phoebus +//======================================================================================== +// © 2022. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract +// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which +// is operated by Triad National Security, LLC for the U.S. +// Department of Energy/National Nuclear Security Administration. All +// rights in the program are reserved by Triad National Security, LLC, +// and the U.S. Department of Energy/National Nuclear Security +// Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, +// distribute copies to the public, perform publicly and display +// publicly, and to permit others to do so. + +#ifndef UTILS_ROBUST_HPP_ +#define UTILS_ROBUST_HPP_ + +#include "basic_types.hpp" +#include +#include + +namespace robust { +using parthenon::Real; + +template +KOKKOS_FORCEINLINE_FUNCTION constexpr auto LARGE() { + return 0.1 * std::numeric_limits::max(); +} +template +KOKKOS_FORCEINLINE_FUNCTION constexpr auto SMALL() { + return 10 * std::numeric_limits::min(); +} +template +KOKKOS_FORCEINLINE_FUNCTION constexpr auto EPS() { + return 10 * std::numeric_limits::epsilon(); +} + +template +KOKKOS_FORCEINLINE_FUNCTION auto make_positive(const T val) { + return std::max(val, EPS()); +} + +KOKKOS_FORCEINLINE_FUNCTION +Real make_bounded(const Real val, const Real vmin, const Real vmax) { + return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); +} + +template +KOKKOS_INLINE_FUNCTION int sgn(const T &val) { + return (T(0) <= val) - (val < T(0)); +} +template +KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) { + const B sgn = b >= 0 ? 1 : -1; + return a / (b + sgn * SMALL()); +} +} // namespace robust + +#endif // UTILS_ROBUST_HPP_ From d2dba04a69940f5b8ffd6fff2eb98e5369ccbae1 Mon Sep 17 00:00:00 2001 From: evan1022 Date: Mon, 11 Mar 2024 13:59:55 -0700 Subject: [PATCH 25/30] Update tracers.cpp fixed a bug in computing the correlations --- src/tracers/tracers.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 65e53225..092c49a0 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -3,7 +3,7 @@ // Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. // Licensed under the BSD 3-Clause License (the "LICENSE"). //======================================================================================== -// Tracer implementation refacored from https://github.com/lanl/phoebus +// Tracer implementation refactored from https://github.com/lanl/phoebus //======================================================================================== // © 2021-2023. Triad National Security, LLC. All rights reserved. // This program was produced under U.S. Government contract @@ -263,7 +263,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // correlations for (s_idx = 0; s_idx < n_lookback; s_idx++) { Kokkos::atomic_add(&corr(0, s_idx), s(0, n) * s(s_idx, n)); - Kokkos::atomic_add(&corr(1, s_idx), sdot(0, n) * s(s_idx, n)); + Kokkos::atomic_add(&corr(1, s_idx), sdot(0, n) * sdot(s_idx, n)); } } }); From 60e8bc1e895052cde1762f666a14c8ef9d94ce44 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 12 Mar 2024 11:29:54 +0100 Subject: [PATCH 26/30] Fix csv header --- src/tracers/tracers.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index b1fff352..0f19fee8 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -309,7 +309,13 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // On startup, write header if (current_cycle == 0) { outfile.open(fname, std::ofstream::out); - outfile << "# Hello world\n"; + outfile << "# cycle, time"; + for (const auto &var : {"corr_s", "corr_sdot", "t_lookback"}) { + for (int i = 0; i < n_lookback; i++) { + outfile << ", " << var << "[" << i << "]"; + } + outfile << std::endl; + } } else { outfile.open(fname, std::ofstream::out | std::ofstream::app); } From 79ddc4f9d5e73421e08ffb3ee6c5675d078f5d04 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 13 Mar 2024 11:51:35 +0100 Subject: [PATCH 27/30] Add s and sdot to csv output --- src/tracers/tracers.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 0f19fee8..88759b32 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -199,7 +199,8 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // TODO(pgrete) Benchmark atomic and potentially update to proper reduction instead of // atomics. // Used for the parallel reduction. Could be reused but this way it's initalized to 0. - parthenon::ParArray2D corr("tracer correlations", 2, n_lookback); + // n_lookback + 1 as it also carries and + parthenon::ParArray2D corr("tracer correlations", 2, n_lookback + 1); int64_t num_particles_total = 0; // Get hydro/mhd fluid vars over all blocks @@ -276,6 +277,8 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { Kokkos::atomic_add(&corr(0, s_idx), s(0, n) * s(s_idx, n)); Kokkos::atomic_add(&corr(1, s_idx), sdot(0, n) * sdot(s_idx, n)); } + Kokkos::atomic_add(&corr(0, n_lookback), s(0, n)); + Kokkos::atomic_add(&corr(1, n_lookback), sdot(0, n)); } }); num_particles_total += swarm->GetNumActive(); @@ -298,7 +301,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { #endif if (parthenon::Globals::my_rank == 0) { // Turn sum into mean - for (int i = 0; i < n_lookback; i++) { + for (int i = 0; i < n_lookback + 1; i++) { corr_h(0, i) /= static_cast(num_particles_total); corr_h(1, i) /= static_cast(num_particles_total); } @@ -309,7 +312,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // On startup, write header if (current_cycle == 0) { outfile.open(fname, std::ofstream::out); - outfile << "# cycle, time"; + outfile << "# cycle, time, s, sdot"; for (const auto &var : {"corr_s", "corr_sdot", "t_lookback"}) { for (int i = 0; i < n_lookback; i++) { outfile << ", " << var << "[" << i << "]"; @@ -322,6 +325,10 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { outfile << tm.ncycle << "," << tm.time; + // and + outfile << "," << corr_h(0, n_lookback); + outfile << "," << corr_h(1, n_lookback); + // and for (int j = 0; j < 2; j++) { for (int i = 0; i < n_lookback; i++) { outfile << "," << corr_h(j, i); From 4e4581a38f8bb3fab0205b899d1f12dd6e46ebbf Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 18 Jul 2024 11:29:29 +0200 Subject: [PATCH 28/30] Fix interface --- src/pgen/turbulence.cpp | 6 +++--- src/tracers/tracers.cpp | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index b825590b..05e2be92 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -361,9 +361,9 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { // Create new particles and get accessor auto new_particles_context = swarm->AddEmptyParticles(num_tracers_per_block); - auto &x = swarm->Get("x").Get(); - auto &y = swarm->Get("y").Get(); - auto &z = swarm->Get("z").Get(); + auto &x = swarm->Get(swarm_position::x::name()).Get(); + auto &y = swarm->Get(swarm_position::y::name()).Get(); + auto &z = swarm->Get(swarm_position::z::name()).Get(); auto &id = swarm->Get("id").Get(); auto swarm_d = swarm->GetDeviceContext(); diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 88759b32..352361bc 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -115,9 +115,9 @@ TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { auto &sd = pmb->swarm_data.Get(); auto &swarm = sd->Get("tracers"); - auto &x = swarm->Get("x").Get(); - auto &y = swarm->Get("y").Get(); - auto &z = swarm->Get("z").Get(); + auto &x = swarm->Get(swarm_position::x::name()).Get(); + auto &y = swarm->Get(swarm_position::y::name()).Get(); + auto &z = swarm->Get(swarm_position::z::name()).Get(); auto &vel_x = swarm->Get("vel_x").Get(); auto &vel_y = swarm->Get("vel_y").Get(); auto &vel_z = swarm->Get("vel_z").Get(); @@ -212,9 +212,9 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { // TODO(pgrete) cleanup once get swarm packs (currently in development upstream) // pull swarm vars - auto &x = swarm->Get("x").Get(); - auto &y = swarm->Get("y").Get(); - auto &z = swarm->Get("z").Get(); + auto &x = swarm->Get(swarm_position::x::name()).Get(); + auto &y = swarm->Get(swarm_position::y::name()).Get(); + auto &z = swarm->Get(swarm_position::z::name()).Get(); auto &vel_x = swarm->Get("vel_x").Get(); auto &vel_y = swarm->Get("vel_y").Get(); auto &vel_z = swarm->Get("vel_z").Get(); From 14976c1134192e977ede416df09b36f6a9bfd6b2 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 9 Sep 2024 21:39:21 +0200 Subject: [PATCH 29/30] Use upstream interfaces --- src/CMakeLists.txt | 2 - src/tracers/tracers.cpp | 10 +- src/utils/interpolation.hpp | 245 ------------------------------------ src/utils/robust.hpp | 66 ---------- 4 files changed, 7 insertions(+), 316 deletions(-) delete mode 100644 src/utils/interpolation.hpp delete mode 100644 src/utils/robust.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 22a3bdac..6d106322 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -23,8 +23,6 @@ add_executable( tracers/tracers.hpp utils/few_modes_ft.cpp utils/few_modes_ft.hpp - utils/interpolation.hpp - utils/robust.hpp ) add_subdirectory(pgen) diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 352361bc..1cd4dee1 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -23,17 +23,21 @@ #include #include -#include "../main.hpp" -#include "../utils/interpolation.hpp" +// Parthenon headers #include "basic_types.hpp" #include "interface/metadata.hpp" #include "kokkos_abstraction.hpp" #include "parthenon_array_generic.hpp" -#include "tracers.hpp" #include "utils/error_checking.hpp" +#include "utils/interpolation.hpp" + +// AthenaPK headers +#include "../main.hpp" +#include "tracers.hpp" namespace Tracers { using namespace parthenon::package::prelude; +namespace LCInterp = parthenon::interpolation::cent::linear; std::shared_ptr Initialize(ParameterInput *pin) { auto tracer_pkg = std::make_shared("tracers"); diff --git a/src/utils/interpolation.hpp b/src/utils/interpolation.hpp deleted file mode 100644 index 6185ffc9..00000000 --- a/src/utils/interpolation.hpp +++ /dev/null @@ -1,245 +0,0 @@ -//======================================================================================== -// AthenaPK - a performance portable block structured AMR astrophysical MHD code. -// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. -// Licensed under the BSD 3-Clause License (the "LICENSE"). -//======================================================================================== -// Interpolation copied/refactored from -// https://github.com/lanl/phoebus and https://github.com/lanl/spiner -//======================================================================================== -// © 2022. Triad National Security, LLC. All rights reserved. This -// program was produced under U.S. Government contract -// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which -// is operated by Triad National Security, LLC for the U.S. -// Department of Energy/National Nuclear Security Administration. All -// rights in the program are reserved by Triad National Security, LLC, -// and the U.S. Department of Energy/National Nuclear Security -// Administration. The Government is granted for itself and others -// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide -// license in this material to reproduce, prepare derivative works, -// distribute copies to the public, perform publicly and display -// publicly, and to permit others to do so. - -#ifndef UTILS_INTERPOLATION_HPP_ -#define UTILS_INTERPOLATION_HPP_ - -// Parthenon includes -#include "utils/error_checking.hpp" -#include -#include -#include - -#include "robust.hpp" - -namespace interpolation { - -using namespace parthenon::package::prelude; -using parthenon::Coordinates_t; - -// From https://github.com/lanl/spiner/blob/main/spiner/regular_grid_1d.hpp -// a poor-man's std::pair -struct weights_t { - Real first, second; - KOKKOS_INLINE_FUNCTION Real &operator[](const int i) { - assert(0 <= i && i <= 1); - return i == 0 ? first : second; - } -}; - -/// Base class for providing interpolation methods on uniformly spaced data. -/// Constructor is provided with spacing, number of support points, and desired -/// shift. GetIndicesAndWeights then updates arrays of indices and weights for -/// calculating the interpolated data. These arrays are of size StencilSize(). -/// Data is forced to zero outside the boundaries. -class Interpolation { - public: - KOKKOS_FUNCTION - Interpolation(const int nSupport, const Real dx, const Real shift) - : nSupport_(nSupport), dx_(dx), shift_(shift), ishift_(std::round(shift)) {} - - KOKKOS_INLINE_FUNCTION - virtual void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const {} - KOKKOS_INLINE_FUNCTION - virtual int StencilSize() const { return 0; } - - static constexpr int maxStencilSize = 2; - - protected: - const int nSupport_; - const Real dx_; - Real shift_; - int ishift_; -}; - -class PiecewiseConstant : public Interpolation { - public: - KOKKOS_FUNCTION - PiecewiseConstant(const int nSupport, const Real dx, const Real shift) - : Interpolation(nSupport, dx, shift) {} - - KOKKOS_INLINE_FUNCTION - void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { - idx[0] = i + ishift_; - wgt[0] = 1.; - if (idx[0] < 0 || idx[0] >= nSupport_) { - idx[0] = 0; - wgt[0] = 0.; - } - } - - KOKKOS_INLINE_FUNCTION - int StencilSize() const override { return 1; } -}; - -class Linear : public Interpolation { - public: - KOKKOS_FUNCTION - Linear(const int nSupport, const Real dx, const Real shift) - : Interpolation(nSupport, dx, shift) { - PARTHENON_FAIL("Not written yet!"); - } - - KOKKOS_INLINE_FUNCTION - void GetIndicesAndWeights(const int i, int *idx, Real *wgt) const override { - idx[0] = std::floor(i + shift_); - idx[1] = idx[0] + 1; - - wgt[0] = wgt[1] = 1. - wgt[0]; - - for (int nsup = 0; nsup < 2; nsup++) { - if (idx[nsup] < 0 || idx[nsup] >= nSupport_) { - idx[nsup] = 0; - wgt[nsup] = 0.; - } - } - } - - KOKKOS_INLINE_FUNCTION - int StencilSize() const override { return 2; } -}; - -// TODO(JMM): Is this interpolation::Do syntax reasonable? An -// alternative path would be a class called "LCInterp with all -// static functions. Then it could have an `operator()` which would -// be maybe nicer? -// TODO(JMM): Merge this w/ what Ben has done. -namespace Cent { -namespace Linear { - -/* - * Get interpolation weights for linear interpolation - * PARAM[IN] - x - location to interpolate to - * PARAM[IN] - nx - number of points along this direction. Used for sanity checks. - * PARAM[IN] - coords - parthenon coords object - * PARAM[OUT] - ix - index of points to interpolate - * PARAM[OUT] - w - weights - */ -template -KOKKOS_INLINE_FUNCTION void GetWeights(const Real x, const int nx, - const Coordinates_t &coords, int &ix, - weights_t &w) { - PARTHENON_DEBUG_REQUIRE( - typeid(Coordinates_t) == typeid(UniformCartesian), - "Interpolation routines currently only work for UniformCartesian"); - const Real min = coords.Xc(0); // assume uniform Cartesian - const Real dx = coords.CellWidthFA(DIR); - ix = std::min(std::max(0, static_cast(robust::ratio(x - min, dx))), nx - 2); - const Real floor = min + ix * dx; - w[1] = robust::ratio(x - floor, dx); - w[0] = 1. - w[1]; -} - -/* - * Trilinear interpolation on a variable or meshblock pack - * PARAM[IN] - b - Meshblock index - * PARAM[IN] - X1, X2, X3 - Coordinate locations - * PARAM[IN] - p - Variable or MeshBlockPack - * PARAM[IN] - v - variable index - */ -template -KOKKOS_INLINE_FUNCTION Real Do3D(int b, const Real X1, const Real X2, const Real X3, - const Pack &p, int v) { - const auto &coords = p.GetCoords(b); - int ix[3]; - weights_t w[3]; - GetWeights(X1, p.GetDim(1), coords, ix[0], w[0]); - GetWeights(X2, p.GetDim(2), coords, ix[1], w[1]); - GetWeights(X3, p.GetDim(3), coords, ix[2], w[2]); - return (w[2][0] * (w[1][0] * (w[0][0] * p(b, v, ix[2], ix[1], ix[0]) + - w[0][1] * p(b, v, ix[2], ix[1], ix[0] + 1)) + - w[1][1] * (w[0][0] * p(b, v, ix[2], ix[1] + 1, ix[0]) + - w[0][1] * p(b, v, ix[2], ix[1] + 1, ix[0] + 1))) + - w[2][1] * (w[1][0] * (w[0][0] * p(b, v, ix[2] + 1, ix[1], ix[0]) + - w[0][1] * p(b, v, ix[2] + 1, ix[1], ix[0] + 1)) + - w[1][1] * (w[0][0] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0]) + - w[0][1] * p(b, v, ix[2] + 1, ix[1] + 1, ix[0] + 1)))); -} - -/* - * Bilinear interpolation on a variable or meshblock pack - * PARAM[IN] - b - Meshblock index - * PARAM[IN] - X1, X2 - Coordinate locations - * PARAM[IN] - p - Variable or MeshBlockPack - * PARAM[IN] - v - variable index - */ -template -KOKKOS_INLINE_FUNCTION Real Do2D(int b, const Real X1, const Real X2, const Pack &p, - int v) { - const auto &coords = p.GetCoords(b); - int ix1, ix2; - weights_t w1, w2; - GetWeights(X1, p.GetDim(1), coords, ix1, w1); - GetWeights(X2, p.GetDim(2), coords, ix2, w2); - return (w2[0] * (w1[0] * p(b, v, 0, ix2, ix1) + w1[1] * p(b, v, 0, ix2, ix1 + 1)) + - w2[1] * - (w1[0] * p(b, v, 0, ix2 + 1, ix1) + w1[1] * p(b, v, 0, ix2 + 1, ix1 + 1))); -} - -/* - * Linear interpolation on a variable or meshblock pack - * PARAM[IN] - b - Meshblock index - * PARAM[IN] - X1 - Coordinate location - * PARAM[IN] - p - Variable or MeshBlockPack - * PARAM[IN] - v - variable index - */ -template -KOKKOS_INLINE_FUNCTION Real Do1D(int b, const Real X1, const Pack &p, int v) { - const auto &coords = p.GetCoords(b); - int ix; - weights_t w; - GetWeights(X1, p.GetDim(1), coords, ix, w); - return w[0] * p(b, v, 0, 0, ix) + w[1] * p(b, v, 0, 0, ix + 1); -} - -/* - * Trilinear or bilinear interpolation on a variable or meshblock pack - * PARAM[IN] - axisymmetric - * PARAM[IN] - b - Meshblock index - * PARAM[IN] - X1, X2, X3 - Coordinate locations - * PARAM[IN] - p - Variable or MeshBlockPack - * PARAM[IN] - v - variable index - */ -// JMM: I know this won't vectorize because of the switch, but it -// probably won't anyway, since we're doing trilinear -// interpolation, which will kill memory locality. Doing it this -// way means we can do trilinear vs bilinear which I think is a -// sufficient win at minimum code bloat. -template -KOKKOS_INLINE_FUNCTION Real Do(int b, const Real X1, const Real X2, const Real X3, - const Pack &p, int v) { - if (p.GetDim(3) > 1) { - return Do3D(b, X1, X2, X3, p, v); - } else if (p.GetDim(2) > 1) { - return Do2D(b, X1, X2, p, v); - } else { // 1D - return Do1D(b, X1, p, v); - } -} - -} // namespace Linear -} // namespace Cent -} // namespace interpolation - -// Convenience Namespace Alias -namespace LCInterp = interpolation::Cent::Linear; - -#endif // UTILS_INTERPOLATION_HPP_ diff --git a/src/utils/robust.hpp b/src/utils/robust.hpp deleted file mode 100644 index 7488abea..00000000 --- a/src/utils/robust.hpp +++ /dev/null @@ -1,66 +0,0 @@ - -//======================================================================================== -// AthenaPK - a performance portable block structured AMR astrophysical MHD code. -// Copyright (c) 2024, Athena-Parthenon Collaboration. All rights reserved. -// Licensed under the BSD 3-Clause License (the "LICENSE"). -//======================================================================================== -// Copied from https://github.com/lanl/phoebus -//======================================================================================== -// © 2022. Triad National Security, LLC. All rights reserved. This -// program was produced under U.S. Government contract -// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which -// is operated by Triad National Security, LLC for the U.S. -// Department of Energy/National Nuclear Security Administration. All -// rights in the program are reserved by Triad National Security, LLC, -// and the U.S. Department of Energy/National Nuclear Security -// Administration. The Government is granted for itself and others -// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide -// license in this material to reproduce, prepare derivative works, -// distribute copies to the public, perform publicly and display -// publicly, and to permit others to do so. - -#ifndef UTILS_ROBUST_HPP_ -#define UTILS_ROBUST_HPP_ - -#include "basic_types.hpp" -#include -#include - -namespace robust { -using parthenon::Real; - -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto LARGE() { - return 0.1 * std::numeric_limits::max(); -} -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto SMALL() { - return 10 * std::numeric_limits::min(); -} -template -KOKKOS_FORCEINLINE_FUNCTION constexpr auto EPS() { - return 10 * std::numeric_limits::epsilon(); -} - -template -KOKKOS_FORCEINLINE_FUNCTION auto make_positive(const T val) { - return std::max(val, EPS()); -} - -KOKKOS_FORCEINLINE_FUNCTION -Real make_bounded(const Real val, const Real vmin, const Real vmax) { - return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); -} - -template -KOKKOS_INLINE_FUNCTION int sgn(const T &val) { - return (T(0) <= val) - (val < T(0)); -} -template -KOKKOS_INLINE_FUNCTION auto ratio(const A &a, const B &b) { - const B sgn = b >= 0 ? 1 : -1; - return a / (b + sgn * SMALL()); -} -} // namespace robust - -#endif // UTILS_ROBUST_HPP_ From da769ad32cf8d99955a3141acec7f6229ff42c7b Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 9 Sep 2024 21:52:44 +0200 Subject: [PATCH 30/30] Fix swarm interface --- src/hydro/hydro_driver.cpp | 6 +++--- src/pgen/turbulence.cpp | 3 +-- src/tracers/tracers.cpp | 10 +++++----- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index a527fb8f..8a9db209 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -15,8 +15,8 @@ #include "basic_types.hpp" #include "bvals/comms/bvals_in_one.hpp" #include "prolong_restrict/prolong_restrict.hpp" -#include #include "utils/error_checking.hpp" +#include // AthenaPK headers #include "../eos/adiabatic_hydro.hpp" #include "../pgen/cluster/agn_triggering.hpp" @@ -705,7 +705,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { { for (auto &pmb : blocks) { auto &tl = sync_region_tr[0]; - auto &sd = pmb->swarm_data.Get(); + auto &sd = pmb->meshblock_data.Get()->GetSwarmData(); auto reset_comms = tl.AddTask(none, &SwarmContainer::ResetCommunication, sd.get()); } @@ -715,7 +715,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { for (int n = 0; n < blocks.size(); n++) { auto &tl = async_region_tr[n]; auto &pmb = blocks[n]; - auto &sd = pmb->swarm_data.Get(); + auto &sd = pmb->meshblock_data.Get()->GetSwarmData(); auto &mbd0 = pmb->meshblock_data.Get("base"); auto tracer_advect = tl.AddTask(none, Tracers::AdvectTracers, mbd0.get(), integrator->dt); diff --git a/src/pgen/turbulence.cpp b/src/pgen/turbulence.cpp index 05e2be92..89269a3a 100644 --- a/src/pgen/turbulence.cpp +++ b/src/pgen/turbulence.cpp @@ -341,8 +341,7 @@ void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData *md) { for (int b = 0; b < md->NumBlocks(); b++) { pmb = md->GetBlockData(b)->GetBlockPointer(); - auto &sd = pmb->swarm_data.Get(); - auto &swarm = pmb->swarm_data.Get()->Get("tracers"); + auto &swarm = pmb->meshblock_data.Get()->GetSwarmData()->Get("tracers"); auto rng_pool = tracer_pkg->Param("rng_pool"); const Real &x_min = pmb->coords.Xf<1>(ib.s); diff --git a/src/tracers/tracers.cpp b/src/tracers/tracers.cpp index 1cd4dee1..42920d4e 100644 --- a/src/tracers/tracers.cpp +++ b/src/tracers/tracers.cpp @@ -84,9 +84,9 @@ std::shared_ptr Initialize(ParameterInput *pin) { // one, but we should check if there's direct way to access Params of other packages. const bool mhd = pin->GetString("hydro", "fluid") == "glmmhd"; - PARTHENON_REQUIRE_THROWS(pin->GetOrAddString("parthenon/mesh", "refinement", "none") == - "none", - "Tracers/swarms currently only supported on uniform meshes."); + PARTHENON_REQUIRE_THROWS( + pin->GetString("parthenon/mesh", "refinement") != "adaptive", + "Tracers/swarms currently only supported on non-adaptive meshes."); if (mhd) { tracer_pkg->AddSwarmValue("B_x", swarm_name, real_swarmvalue_metadata); @@ -116,7 +116,7 @@ std::shared_ptr Initialize(ParameterInput *pin) { TaskStatus AdvectTracers(MeshBlockData *mbd, const Real dt) { auto *pmb = mbd->GetParentPointer(); - auto &sd = pmb->swarm_data.Get(); + auto &sd = pmb->meshblock_data.Get()->GetSwarmData(); auto &swarm = sd->Get("tracers"); auto &x = swarm->Get(swarm_position::x::name()).Get(); @@ -211,7 +211,7 @@ TaskStatus FillTracers(MeshData *md, parthenon::SimTime &tm) { const auto &prim_pack = md->PackVariables(std::vector{"prim"}); for (int b = 0; b < md->NumBlocks(); b++) { auto *pmb = md->GetBlockData(b)->GetBlockPointer(); - auto &sd = pmb->swarm_data.Get(); + auto &sd = pmb->meshblock_data.Get()->GetSwarmData(); auto &swarm = sd->Get("tracers"); // TODO(pgrete) cleanup once get swarm packs (currently in development upstream)