From 8720309cefe34d3488cd4c620d6f73cb734d40b5 Mon Sep 17 00:00:00 2001 From: Shengtai Li Date: Fri, 17 Oct 2025 13:52:55 -0600 Subject: [PATCH 1/3] commit change for user-src --- src/artemis.hpp | 2 ++ src/artemis_driver.cpp | 5 ++- src/pgen/disk.hpp | 68 +++++++++++++++++++++++++++++++++++ src/pgen/problem_modifier.hpp | 4 ++- src/utils/artemis_utils.hpp | 10 ++++++ 5 files changed, 87 insertions(+), 2 deletions(-) diff --git a/src/artemis.hpp b/src/artemis.hpp index 89e61c68..de4830fe 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -238,6 +238,8 @@ inline int ProblemDimension(parthenon::ParameterInput *pin) { // Custom AMR criteria namespace artemis { extern std::function *mbd)> ProblemCheckRefinementBlock; +extern std::function *md, const Real time, const Real dt)> + UserSourceTerm; } // namespace artemis #endif // ARTEMIS_ARTEMIS_HPP_ diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 8efba20e..c9353908 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -302,9 +302,12 @@ TaskCollection ArtemisDriver::StepTasks() { tl.AddTask(rt_src, RotatingFrame::RotatingFrameForce, u0.get(), time, bdt); } + auto user_src = + tl.AddTask(rframe_src, ArtemisUtils::ProblemSourceTerm, 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); } diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index b9eb307f..30baf8b8 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -890,6 +890,74 @@ void DiskBoundaryExtrap(std::shared_ptr> &mbd, bool coarse) }); } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus ProblemGeneratorSourceTerm() +//! \brief Custom source term for disk pgen +template +TaskStatus UserSourceTerm(MeshData *md, const Real time, const Real dt) { + if (GEOM == Coordinates::spherical3D || GEOM == Coordinates::spherical2D) { + // reset the gas and dust state above nH_init to initial values + 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"); + + static auto desc = + MakePackDescriptor(resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + // Disk parameters + const auto &pgen = artemis_pkg->template Param("disk_params"); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "UniformGravity", 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) { + // Extract coordinates + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const auto &hx = coords.GetScaleFactors(); + const auto &xv = coords.GetCellCenter(); + const auto &xcyl = coords.ConvertToCyl(xv); + const Real H = xcyl[0] * pgen.h0 * std::pow(xcyl[0] / pgen.r0, pgen.flare); + if (std::abs(xcyl[2]) > 5.0 * H) { + if (do_gas) { + for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) = 0.0; + vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) = 0.0; + vmesh(b, gas::cons::total_energy(n), k, j, i) -= + 0.5 * + (SQR(vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i)) + + SQR(vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i))) / + vmesh(b, gas::cons::density(n), k, j, i); + } + } + + if (do_dust) { + for (int n = 0; n < vmesh.GetSize(b, dust::prim::density()); ++n) { + vmesh(b, dust::cons::density(n), k, j, i) = 0.0; + vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) = 0.0; + vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) = 0.0; + vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) = 0.0; + } + } + } // end if (std::abs(xcyl[2]) > 5.0 * H) + }); + } + + return TaskStatus::complete; +} + //---------------------------------------------------------------------------------------- //! \fn AmrTag ProblemCheckRefinementBlock() //! \brief Refinement criterion for disk pgen diff --git a/src/pgen/problem_modifier.hpp b/src/pgen/problem_modifier.hpp index 67f3075d..c0de64f0 100644 --- a/src/pgen/problem_modifier.hpp +++ b/src/pgen/problem_modifier.hpp @@ -30,7 +30,8 @@ namespace artemis { std::function *mbd)> ProblemCheckRefinementBlock = nullptr; - +std::function *md, const Real time, const Real dt)> + UserSourceTerm = nullptr; } // namespace artemis // Problem modifiers @@ -66,6 +67,7 @@ void ProblemModifier(parthenon::ParthenonManager *pman) { } else if (artemis_problem == "disk") { artemis::ProblemCheckRefinementBlock = disk::ProblemCheckRefinementBlock; + artemis::UserSourceTerm = disk::UserSourceTerm; pman->app_input->RegisterBoundaryCondition(BF::inner_x1, "ic", disk::DiskBoundaryIC); diff --git a/src/utils/artemis_utils.hpp b/src/utils/artemis_utils.hpp index de211f36..2ad87f14 100644 --- a/src/utils/artemis_utils.hpp +++ b/src/utils/artemis_utils.hpp @@ -35,6 +35,16 @@ KOKKOS_FORCEINLINE_FUNCTION Real VDot(const V1 &a, const V2 &b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } +//---------------------------------------------------------------------------------------- +//! \fn int ArtemisUtils::ProblemSourceTerm +//! \brief Wrapper function for user-defined source term, checking for nullptr +static TaskStatus ProblemSourceTerm(MeshData *md, const Real time, const Real dt) { + if (artemis::UserSourceTerm != nullptr) { + return artemis::UserSourceTerm(md, time, dt); + } + return TaskStatus::complete; +} + //---------------------------------------------------------------------------------------- //! \fn Real ArtemisUtils::DualEnergySIE(vmesh, const int b, const int n, const int k, //! const int j, const int i, const Real de_switch, From 5640aaaa3ed2fe5efc506a489ada63899a391831 Mon Sep 17 00:00:00 2001 From: Shengtai Li Date: Tue, 13 Jan 2026 11:53:48 -0700 Subject: [PATCH 2/3] change names to be consistent --- src/pgen/disk.hpp | 49 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index 30baf8b8..e4b31ea5 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -1,3 +1,4 @@ + //======================================================================================== // (C) (or copyright) 2023-2025. Triad National Security, LLC. All rights reserved. // @@ -65,6 +66,7 @@ struct DiskParams { Real mdot; Real temp_soft2; Real kbmu, ar; + Real nH_reset; bool do_gas, do_dust, do_moment, do_imc; bool nbody_temp; bool quiet_start; @@ -257,6 +259,7 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) { disk_params.l0 = pin->GetOrAddReal("problem", "l0", 0.0); disk_params.dust_to_gas = pin->GetOrAddReal("problem", "dust_to_gas", 0.01); disk_params.temp_soft2 = pin->GetOrAddReal("problem", "temp_soft", 0.0); + disk_params.nH_reset = pin->GetOrAddReal("problem", "nH_reset", 0.0); const auto mu = gas_pkg->Param("mu"); auto &constants = artemis_pkg->Param("constants"); const auto &eos = gas_pkg->Param("eos_h"); @@ -891,17 +894,27 @@ void DiskBoundaryExtrap(std::shared_ptr> &mbd, bool coarse) } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus ProblemGeneratorSourceTerm() +//! \fn TaskStatus UserSourceTerm() //! \brief Custom source term for disk pgen template TaskStatus UserSourceTerm(MeshData *md, const Real time, const Real dt) { if (GEOM == Coordinates::spherical3D || GEOM == Coordinates::spherical2D) { - // reset the gas and dust state above nH_init to initial values + // reset the gas and dust state above nH_reset: gas to initial values, dust to zeros using parthenon::MakePackDescriptor; auto pm = md->GetParentPointer(); - auto &resolved_pkgs = pm->resolved_packages; auto &artemis_pkg = pm->packages.Get("artemis"); + + // Disk parameters + const auto &pgen = artemis_pkg->template Param("disk_params"); + if (pgen.nH_reset <= 0.0) return TaskStatus::complete; + + // Extract gas package and params + auto &gas_pkg = pm->packages.Get("gas"); + const auto &eos_d = gas_pkg->template Param("eos_d"); + + auto &resolved_pkgs = pm->resolved_packages; + const bool do_gas = artemis_pkg->template Param("do_gas"); const bool do_dust = artemis_pkg->template Param("do_dust"); @@ -917,29 +930,35 @@ TaskStatus UserSourceTerm(MeshData *md, const Real time, const Real dt) { const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); - // Disk parameters - const auto &pgen = artemis_pkg->template Param("disk_params"); - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "UniformGravity", parthenon::DevExecSpace(), 0, + DEFAULT_LOOP_PATTERN, "UserSourceTerm", 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) { // Extract coordinates geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); const auto &hx = coords.GetScaleFactors(); const auto &xv = coords.GetCellCenter(); - const auto &xcyl = coords.ConvertToCyl(xv); + const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); const Real H = xcyl[0] * pgen.h0 * std::pow(xcyl[0] / pgen.r0, pgen.flare); - if (std::abs(xcyl[2]) > 5.0 * H) { + if (std::abs(xcyl[2]) > pgen.nH_reset * H) { if (do_gas) { + const Real nu = ViscosityProfile(pgen, eos_d, xcyl[0], xcyl[2]); + const Real vr = pgen.quiet_start ? 0.0 : -1.5 * nu / xcyl[0]; for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { - vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) = 0.0; - vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) = 0.0; - vmesh(b, gas::cons::total_energy(n), k, j, i) -= + const Real etot_res = + vmesh(b, gas::cons::total_energy(n), k, j, i) - + 0.5 * + (SQR(vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i)) + + SQR(vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i))) / + vmesh(b, gas::cons::density(n), k, j, i); + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) = vr * xcyl[0] / xv[0]; + vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) = vr * xcyl[2] / xv[0]; + vmesh(b, gas::cons::total_energy(n), k, j, i) = + etot_res + 0.5 * - (SQR(vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i)) + - SQR(vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i))) / - vmesh(b, gas::cons::density(n), k, j, i); + (SQR(vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i)) + + SQR(vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i))) / + vmesh(b, gas::cons::density(n), k, j, i); } } From 16c5d84c3b7ff7a49ca58963409fd53d76a25485 Mon Sep 17 00:00:00 2001 From: Shengtai Li Date: Tue, 13 Jan 2026 12:48:07 -0700 Subject: [PATCH 3/3] add a parameter in formula of stopping times for Stokes regime --- src/drag/drag.hpp | 14 ++++++++++---- src/pgen/disk.hpp | 8 ++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/drag/drag.hpp b/src/drag/drag.hpp index 4358cdd8..ae9fd1b4 100644 --- a/src/drag/drag.hpp +++ b/src/drag/drag.hpp @@ -125,6 +125,7 @@ struct StoppingTimeParams { DragModel model; ParArray1D tau; Real tau_max, tau_min; + Real x1_power; StoppingTimeParams(std::string block_name, ParameterInput *pin) { const std::string choice = pin->GetString(block_name, "type"); @@ -146,6 +147,7 @@ struct StoppingTimeParams { scale = pin->GetOrAddReal(block_name, "scale", 1.0); tau_max = pin->GetOrAddReal(block_name, "maximum", 1e99); tau_min = pin->GetOrAddReal(block_name, "minimum", 0.0); + x1_power = pin->GetOrAddReal(block_name, "x1_power", 0.0); auto h_tau = tau.GetHostMirror(); for (int n = 0; n < nd; n++) { h_tau(n) = scale; @@ -463,9 +465,11 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real // Extract Stokes specific parameters [[maybe_unused]] auto &grain_density_ = grain_density; [[maybe_unused]] Real vth = Null(); + Real x1fac = 1.0; if constexpr (DRAG == DragModel::stokes) { const Real gm1 = eos_d.GruneisenParamFromDensityInternalEnergy(dg, sieg); vth = std::sqrt(8.0 / M_PI * gm1 * sieg); + x1fac = std::pow(xv[0], tp.x1_power); } // First pass to collect \sum rho' and \sum rho' v and compute new vg @@ -493,8 +497,9 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real Real tc = tp.tau(id); [[maybe_unused]] auto &sizes_ = sizes; if constexpr (DRAG == DragModel::stokes) { - tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / - dg * sizes_(id) / vth)); + tc = + std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / dg * + sizes_(id) / vth * x1fac)); } const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); for (int d = 0; d < 3; d++) { @@ -529,8 +534,9 @@ TaskStatus SimpleDragSourceImpl(MeshData *md, const Real time, const Real Real tc = tp.tau(id); [[maybe_unused]] auto &sizes_ = sizes; if constexpr (DRAG == DragModel::stokes) { - tc = std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / - dg * sizes_(id) / vth)); + tc = + std::max(tp.tau_min, std::min(tp.tau_max, tp.scale * grain_density_ / dg * + sizes_(id) / vth * x1fac)); } const Real alpha = dt * ((tc <= 0.0) ? Big() : 1.0 / tc); // Update dust momenta diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp index e4b31ea5..a28da111 100644 --- a/src/pgen/disk.hpp +++ b/src/pgen/disk.hpp @@ -203,7 +203,11 @@ KOKKOS_INLINE_FUNCTION State ComputeDiskProfile( const Real vp2 = vk2 + (dpdr / res.gdens) * xcyl[0]; const Real vp = (vp2 < 0.0) ? 0.0 : std::sqrt(vp2); const Real nu = ViscosityProfile(pgen, eos_d, rt, xcyl[2]); - const Real vr = pgen.quiet_start ? 0.0 : -1.5 * nu / xcyl[0]; + const Real vr_visc = pgen.alpha > 0.0 + ? (-3.0 * pgen.alpha / std::sqrt(xcyl[0]) / res.gdens * + (SQR(xcyl[0]) * dpdr + pres * 2.0 * xcyl[0])) + : (-1.5 * nu / xcyl[0]); + const Real vr = pgen.quiet_start ? 0.0 : vr_visc; // Construct the total cylindrical velocity const Real vcyl[3] = {vr, vp - pgen.omf * xcyl[0], 0.0}; @@ -938,7 +942,7 @@ TaskStatus UserSourceTerm(MeshData *md, const Real time, const Real dt) { geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); const auto &hx = coords.GetScaleFactors(); const auto &xv = coords.GetCellCenter(); - const auto &[xcyl, ex1, ex2, ex3] = coords.ConvertToCylWithVec(xv); + const auto &xcyl = coords.ConvertToCyl(xv); const Real H = xcyl[0] * pgen.h0 * std::pow(xcyl[0] / pgen.r0, pgen.flare); if (std::abs(xcyl[2]) > pgen.nH_reset * H) { if (do_gas) {