From 85ecde79ac53974400c8b858ceda6e65fb3e5123 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Thu, 23 Jul 2026 09:06:28 -0600 Subject: [PATCH 1/2] Add a registry for user source term tasks. Add normal wave-killing to disk pgen as an example. This could be expanded to other registries for tasks in different parts of the operator split in the future --- src/artemis.cpp | 31 ++++++ src/artemis.hpp | 16 +++ src/artemis_driver.cpp | 12 ++- src/pgen/disk.hpp | 194 ++++++++++++++++++++++++++++++++++ src/pgen/problem_modifier.hpp | 7 ++ 5 files changed, 258 insertions(+), 2 deletions(-) diff --git a/src/artemis.cpp b/src/artemis.cpp index 7139c290..9d5f37ce 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -11,6 +11,11 @@ // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== +// This file was created in part by generative AI + +// C++ includes +#include + // Artemis includes #include "artemis.hpp" #include "artemis_driver.hpp" @@ -36,6 +41,32 @@ namespace artemis { +namespace { +std::vector user_source_tasks; +} // namespace + +//---------------------------------------------------------------------------------------- +//! \fn void Artemis::RegisterUserSourceTask +//! \brief Registers a named problem-generator source task in execution order +void RegisterUserSourceTask(const std::string &name, UserSourceTaskFn function) { + PARTHENON_REQUIRE(!name.empty(), "User source task names cannot be empty"); + PARTHENON_REQUIRE(function != nullptr, "Cannot register a null user source task"); + for (const auto &task : user_source_tasks) { + PARTHENON_REQUIRE(task.name != name, "Duplicate user source task name: " + name); + } + user_source_tasks.push_back({name, std::move(function)}); +} + +//---------------------------------------------------------------------------------------- +//! \fn void Artemis::ClearUserSourceTasks +//! \brief Clears source tasks registered by a previous problem modifier +void ClearUserSourceTasks() { user_source_tasks.clear(); } + +//---------------------------------------------------------------------------------------- +//! \fn const std::vector &Artemis::GetUserSourceTasks +//! \brief Returns problem-generator source tasks in registration order +const std::vector &GetUserSourceTasks() { return user_source_tasks; } + //---------------------------------------------------------------------------------------- //! \fn Packages_t Artemis::ProcessPackages //! \brief Process and initialize relevant packages diff --git a/src/artemis.hpp b/src/artemis.hpp index 20e13e93..d59f7dfb 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -10,12 +10,16 @@ // 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. //======================================================================================== + +// This file was created in part by generative AI #ifndef ARTEMIS_ARTEMIS_HPP_ #define ARTEMIS_ARTEMIS_HPP_ // C++ includes +#include #include #include +#include // Parthenon includes #include @@ -235,6 +239,18 @@ inline int ProblemDimension(parthenon::ParameterInput *pin) { // Custom AMR criteria namespace artemis { extern std::function *mbd)> ProblemCheckRefinementBlock; + +using UserSourceTaskFn = + std::function *md, const Real time, const Real dt)>; + +struct UserSourceTask { + std::string name; + UserSourceTaskFn function; +}; + +void RegisterUserSourceTask(const std::string &name, UserSourceTaskFn function); +void ClearUserSourceTasks(); +const std::vector &GetUserSourceTasks(); } // namespace artemis #endif // ARTEMIS_ARTEMIS_HPP_ diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a4af03db..0ab26bb8 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -11,6 +11,8 @@ // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== +// This file was created in part by generative AI + // NOTE(PDM): The following is largely borrowed from the open-source LANL phoebus // software, with additional extensions motivated by other downstream development. @@ -324,11 +326,17 @@ TaskCollection ArtemisDriver::StepTasks() { tl.AddTask(rt_src, RotatingFrame::RotatingFrameForce, u0.get(), time, bdt); } + // Apply problem-generator source terms in registration order + TaskID user_src = rframe_src; + for (const auto &task : GetUserSourceTasks()) { + user_src = tl.AddTask(user_src, task.name, task.function, u0.get(), time, bdt); + } + // Apply drag source term // NOTE(@pdmullen): RK integrated, operator split drag (RHS computed from U) - TaskID drag_src = rframe_src; + TaskID drag_src = user_src; if (do_drag) { - drag_src = tl.AddTask(rframe_src, Drag::DragSource, u0.get(), time, bdt); + drag_src = tl.AddTask(user_src, Drag::DragSource, u0.get(), time, bdt); } // Apply cooling source term diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index fd089728..61f0061a 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -14,6 +14,8 @@ // Copyright(C) 2020 James M. Stone and the Athena code team // Licensed under the 3-clause BSD License (the "LICENSE") //======================================================================================== + +// This file was created in part by generative AI #ifndef PGEN_DISK_HPP_ #define PGEN_DISK_HPP_ //! \file disk.hpp @@ -74,6 +76,30 @@ struct DiskParams { bool multi_d, three_d; }; +//! \struct WaveKillingParams +//! \brief Coordinate-space bounds and rates for disk wave-killing zones +struct WaveKillingParams { + std::array xmin, xmax; + std::array inner, outer; + std::array inner_rate, outer_rate; +}; + +KOKKOS_INLINE_FUNCTION +Real WaveKillingRate(const WaveKillingParams &wave, const std::array &xv) { + Real rate = 0.0; + for (int d = 0; d < 3; ++d) { + if (wave.inner_rate[d] > 0.0 && xv[d] < wave.inner[d]) { + const Real ramp = (xv[d] - wave.inner[d]) / (wave.inner[d] - wave.xmin[d]); + rate += wave.inner_rate[d] * SQR(ramp); + } + if (wave.outer_rate[d] > 0.0 && xv[d] > wave.outer[d]) { + const Real ramp = (xv[d] - wave.outer[d]) / (wave.xmax[d] - wave.outer[d]); + rate += wave.outer_rate[d] * SQR(ramp); + } + } + return rate; +} + struct State { Real gdens = Null(); Real gtemp = Null(); @@ -347,6 +373,39 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { disk_params.nbody_temp = pin->GetOrAddBoolean("problem", "nbody_temp", false); disk_params.nbody_temp = disk_params.nbody_temp && params.Get("do_nbody"); params.Add("disk_params", disk_params); + + if (pin->DoesBlockExist("problem/wave_killing")) { + constexpr char block[] = "problem/wave_killing"; + WaveKillingParams wave; + wave.xmin = {params.Get("x1min"), params.Get("x2min"), + params.Get("x3min")}; + wave.xmax = {params.Get("x1max"), params.Get("x2max"), + params.Get("x3max")}; + + const std::array axis = {"x1", "x2", "x3"}; + for (int d = 0; d < 3; ++d) { + wave.inner[d] = pin->GetOrAddReal(block, axis[d] + "_inner", wave.xmin[d]); + wave.outer[d] = pin->GetOrAddReal(block, axis[d] + "_outer", wave.xmax[d]); + wave.inner_rate[d] = pin->GetOrAddReal(block, axis[d] + "_inner_rate", 0.0); + wave.outer_rate[d] = pin->GetOrAddReal(block, axis[d] + "_outer_rate", 0.0); + + PARTHENON_REQUIRE(wave.inner_rate[d] >= 0.0 && wave.outer_rate[d] >= 0.0, + "Wave-killing rates must be nonnegative"); + PARTHENON_REQUIRE(wave.inner[d] >= wave.xmin[d] && + wave.outer[d] <= wave.xmax[d] && + wave.inner[d] <= wave.outer[d], + "Wave-killing bounds must satisfy xmin <= inner <= outer <= " + "xmax"); + PARTHENON_REQUIRE(wave.inner_rate[d] == 0.0 || wave.inner[d] > wave.xmin[d], + "An active inner wave-killing zone must extend beyond xmin"); + PARTHENON_REQUIRE(wave.outer_rate[d] == 0.0 || wave.outer[d] < wave.xmax[d], + "An active outer wave-killing zone must begin before xmax"); + PARTHENON_REQUIRE(nx[d] > 1 || + (wave.inner_rate[d] == 0.0 && wave.outer_rate[d] == 0.0), + "Wave killing cannot be active in a collapsed dimension"); + } + params.Add("wave_killing_params", wave); + } } } @@ -448,6 +507,141 @@ inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { if (dp.do_imc) jaybenne::InitializeRadiation(md.get(), true); } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus disk::WaveKilling +//! \brief Relaxes disk gas and dust toward their initial profiles near mesh boundaries +template +TaskStatus WaveKilling(MeshData *md, const Real /* time */, const Real dt) { + PARTHENON_INSTRUMENT + using parthenon::MakePackDescriptor; + + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + auto &artemis_pkg = pm->packages.Get("artemis"); + const bool do_gas = artemis_pkg->template Param("do_gas"); + const bool do_dust = artemis_pkg->template Param("do_dust"); + const auto &cpars = artemis_pkg->template Param("coord_params"); + const auto disk_params = artemis_pkg->template Param("disk_params"); + const auto wave = artemis_pkg->template Param("wave_killing_params"); + + auto &gas_pkg = pm->packages.Get("gas"); + const auto &eos_d = gas_pkg->template Param("eos_d"); + const Real de_switch = gas_pkg->template Param("de_switch"); + const Real dflr_gas = gas_pkg->template Param("dfloor"); + const Real sieflr_gas = gas_pkg->template Param("siefloor"); + + Real dflr_dust = Null(); + if (do_dust) { + dflr_dust = pm->packages.Get("dust")->template Param("dfloor"); + } + + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + ParArray1D particles; + int npart = 0; + if (disk_params.nbody_temp) { + auto &nbody_pkg = pm->packages.Get("nbody"); + particles = nbody_pkg->template Param>("particles"); + npart = static_cast(particles.size()); + } + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "DiskWaveKilling", parthenon::DevExecSpace(), 0, + md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &xv = coords.GetCellCenter(vg, b, k, j, i); + const Real rate = WaveKillingRate(wave, xv); + if (rate <= 0.0) return; + + const auto &dx = coords.GetCellWidths(vg, b, k, j, i); + const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); + const Real fac = rate * dt / (1.0 + rate * dt); + const auto target = + ComputeDiskProfile(disk_params, coords, xv, dx, k, j, i, eos_d, do_gas, + do_dust, particles, npart); + + if (do_gas) { + const Real target_vel[3] = {target.gvel1, target.gvel2, target.gvel3}; + for (int n = 0; n < vmesh.GetSize(b, gas::cons::density()); ++n) { + Real &dens = vmesh(b, gas::cons::density(n), k, j, i); + Real &mom1 = vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i); + Real &mom2 = vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i); + Real &mom3 = vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i); + Real &etot = vmesh(b, gas::cons::total_energy(n), k, j, i); + Real &eint = vmesh(b, gas::cons::internal_energy(n), k, j, i); + + dens = std::max(dens, dflr_gas); + Real sie = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx); + sie = std::max(sie, sieflr_gas); + + const Real vel[3] = {mom1 / (dens * hx[0]), mom2 / (dens * hx[1]), + mom3 / (dens * hx[2])}; + const Real temp = eos_d.TemperatureFromDensityInternalEnergy(dens, sie); + const Real new_dens = std::max(dflr_gas, dens + fac * (target.gdens - dens)); + const Real new_temp = temp + fac * (target.gtemp - temp); + const Real new_sie = + std::max(sieflr_gas, + eos_d.InternalEnergyFromDensityTemperature(new_dens, new_temp)); + const Real new_vel[3] = {vel[0] + fac * (target_vel[0] - vel[0]), + vel[1] + fac * (target_vel[1] - vel[1]), + vel[2] + fac * (target_vel[2] - vel[2])}; + + const Real old_ke = 0.5 * dens * (SQR(vel[0]) + SQR(vel[1]) + SQR(vel[2])); + const Real new_ke = + 0.5 * new_dens * (SQR(new_vel[0]) + SQR(new_vel[1]) + SQR(new_vel[2])); + const Real old_eint = dens * sie; + const Real new_eint = new_dens * new_sie; + + dens = new_dens; + mom1 = new_dens * new_vel[0] * hx[0]; + mom2 = new_dens * new_vel[1] * hx[1]; + mom3 = new_dens * new_vel[2] * hx[2]; + eint = new_eint; + etot += (new_eint - old_eint) + (new_ke - old_ke); + etot = std::max(etot, new_eint + new_ke); + } + } + + if (do_dust) { + const Real target_vel[3] = {target.dvel1, target.dvel2, target.dvel3}; + for (int n = 0; n < vmesh.GetSize(b, dust::cons::density()); ++n) { + Real &dens = vmesh(b, dust::cons::density(n), k, j, i); + Real &mom1 = vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i); + Real &mom2 = vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i); + Real &mom3 = vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i); + + dens = std::max(dens, dflr_dust); + const Real vel[3] = {mom1 / (dens * hx[0]), mom2 / (dens * hx[1]), + mom3 / (dens * hx[2])}; + const Real new_dens = std::max(dflr_dust, dens + fac * (target.ddens - dens)); + const Real new_vel[3] = {vel[0] + fac * (target_vel[0] - vel[0]), + vel[1] + fac * (target_vel[1] - vel[1]), + vel[2] + fac * (target_vel[2] - vel[2])}; + + dens = new_dens; + mom1 = new_dens * new_vel[0] * hx[0]; + mom2 = new_dens * new_vel[1] * hx[1]; + mom3 = new_dens * new_vel[2] * hx[2]; + } + } + }); + + return TaskStatus::complete; +} + //---------------------------------------------------------------------------------------- //! \fn void Disk::DiskBoundaryVisc() //! \brief Sets inner or outer X1 boundary condition to the initial condition diff --git a/src/pgen/problem_modifier.hpp b/src/pgen/problem_modifier.hpp index 162e3c02..45ec5961 100644 --- a/src/pgen/problem_modifier.hpp +++ b/src/pgen/problem_modifier.hpp @@ -10,6 +10,8 @@ // 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. //======================================================================================== + +// This file was created in part by generative AI #ifndef PGEN_PROBLEM_MODIFIER_HPP_ #define PGEN_PROBLEM_MODIFIER_HPP_ @@ -44,6 +46,8 @@ void ProblemModifier(parthenon::ParthenonManager *pman) { using BF = parthenon::BoundaryFace; using ID = parthenon::IndexDomain; + ClearUserSourceTasks(); + std::string artemis_problem = pman->pinput->GetOrAddString("artemis", "problem", "unset"); @@ -66,6 +70,9 @@ void ProblemModifier(parthenon::ParthenonManager *pman) { } else if (artemis_problem == "disk") { artemis::ProblemCheckRefinementBlock = disk::ProblemCheckRefinementBlock; + if (pman->pinput->DoesBlockExist("problem/wave_killing")) { + RegisterUserSourceTask("disk_wave_killing", disk::WaveKilling); + } pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "ic", disk::DiskBoundaryIC); From f5354a630ec63b2f8bc51d745ac40ba91b61b0c2 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Thu, 23 Jul 2026 09:08:38 -0600 Subject: [PATCH 2/2] Add wave-killing to disk test --- inputs/disk/wave_killing.in | 107 ++++++++++++++++++++++++++++++++++++ tst/scripts/disk/disk.py | 28 ++++++++++ 2 files changed, 135 insertions(+) create mode 100644 inputs/disk/wave_killing.in diff --git a/inputs/disk/wave_killing.in b/inputs/disk/wave_killing.in new file mode 100644 index 00000000..c70e2058 --- /dev/null +++ b/inputs/disk/wave_killing.in @@ -0,0 +1,107 @@ +# ======================================================================================== +# (C) (or copyright) 2026. 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. +# ======================================================================================== + +# This file was created in part by generative AI +# + +problem = disk +coordinates = axisymmetric + + +problem_id = disk_wave_killing + + +variables = gas.prim.density, & + gas.prim.velocity, & + gas.prim.pressure, & + dust.prim.density, & + dust.prim.velocity +file_type = hdf5 +dt = 1.0 + + +nlim = 3 +tlim = 1.0 +integrator = rk2 +ncycle_out = 1 + + +nx1 = 32 +x1min = 0.5 +x1max = 2.0 +ix1_bc = ic +ox1_bc = ic + +nx2 = 16 +x2min = -0.2 +x2max = 0.2 +ix2_bc = ic +ox2_bc = ic + +nx3 = 1 +x3min = -0.5 +x3max = 0.5 +ix3_bc = periodic +ox3_bc = periodic + + +nx1 = 8 +nx2 = 8 +nx3 = 1 + + +gas = true +dust = true +gravity = true + + +cfl = 0.3 +reconstruct = plm +riemann = hllc-general +dfloor = 1.0e-10 +siefloor = 1.0e-10 + + +gamma = 1.4 + + +cfl = 0.3 +nspecies = 2 +size_input = direct +sizes = 1.0e-4, 1.0e-3 +grain_density = 1.0 +reconstruct = plm +riemann = hlle +dfloor = 1.0e-10 + + +mass = 1.0 + + +r0 = 1.0 +rho0 = 1.0 +dslope = -1.0 +flare = 0.0 +h0 = 0.1 +polytropic_index = 1.0 +quiet_start = true + + +x1_inner = 0.7 +x1_outer = 1.8 +x1_inner_rate = 5.0 +x1_outer_rate = 5.0 +x2_inner = -0.15 +x2_outer = 0.15 +x2_inner_rate = 5.0 +x2_outer_rate = 5.0 diff --git a/tst/scripts/disk/disk.py b/tst/scripts/disk/disk.py index 0d12443d..8b9c91d1 100644 --- a/tst/scripts/disk/disk.py +++ b/tst/scripts/disk/disk.py @@ -95,6 +95,12 @@ def run(**kwargs): ), ) + artemis.run( + _nranks, + "disk/wave_killing.in", + ["parthenon/job/problem_id={}_wave_killing".format(_file_id)], + ) + # Analyze outputs def analyze(): @@ -185,6 +191,28 @@ def analyze(): ) bad |= mybad + wave_file = os.path.join( + artemis.get_data_dir(), "{}_wave_killing.out1.final.phdf".format(_file_id) + ) + with h5py.File(wave_file, "r") as f: + gas_dens = f["gas.prim.density_0"][...] + gas_temp = f["gas.prim.pressure_0"][...] / gas_dens + gas_vel = f["gas.prim.velocity_0"][...] + dust_dens = np.stack( + [f["dust.prim.density_{:d}".format(n)][...] for n in range(2)] + ) + dust_vel = np.stack( + [f["dust.prim.velocity_{:d}".format(n)][...] for n in range(2)] + ) + + wave_fields = (gas_dens, gas_temp, gas_vel, dust_dens, dust_vel) + wave_bad = any(np.any(~np.isfinite(field)) for field in wave_fields) + wave_bad |= np.any(gas_dens <= 0.0) or np.any(gas_temp <= 0.0) + wave_bad |= np.any(dust_dens <= 0.0) + if wave_bad: + logger.debug("disk_wave_killing FAILED: non-finite or non-positive state") + bad |= wave_bad + return not bad