From 2a973d6adb3112a628d5a037afb29bdb6252e1f0 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 16 Dec 2025 16:17:57 -0500 Subject: [PATCH 1/8] Debugging self-gravity --- inputs/self-gravity/poly.in | 103 +++++++++ src/CMakeLists.txt | 6 +- src/artemis.cpp | 1 + src/artemis.hpp | 5 + src/artemis_driver.cpp | 7 +- src/gravity/gravity.cpp | 199 ++++++++++++++++- src/gravity/gravity.hpp | 13 +- src/gravity/poisson_driver.cpp | 77 +++++++ src/gravity/poisson_equation.hpp | 222 +++++++++++++++++++ src/gravity/self_gravity.cpp | 145 ++++++++++++ src/pgen/pgen.hpp | 3 + src/pgen/polytrope.hpp | 113 ++++++++++ src/utils/integrators/artemis_integrator.hpp | 26 +-- 13 files changed, 900 insertions(+), 20 deletions(-) create mode 100644 inputs/self-gravity/poly.in create mode 100644 src/gravity/poisson_driver.cpp create mode 100644 src/gravity/poisson_equation.hpp create mode 100644 src/gravity/self_gravity.cpp create mode 100644 src/pgen/polytrope.hpp diff --git a/inputs/self-gravity/poly.in b/inputs/self-gravity/poly.in new file mode 100644 index 00000000..75ebf774 --- /dev/null +++ b/inputs/self-gravity/poly.in @@ -0,0 +1,103 @@ +# ======================================================================================== +# (C) (or copyright) 2023-2024. 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. +# ======================================================================================== + + +problem = polytrope # name of the pgen +coordinates = cartesian # coordinate system + + +problem_id = polytrope # problem ID: basename of output filenames + + +file_type = hdf5 # HDF5 data dump +variables = gas.prim.density, & + gas.prim.velocity, & + gas.prim.pressure, & + grav.phi, & + grav.rhs +dt = 1.0e-100 # time increment between outputs + + +file_type = hst # HDF5 data dump +dt = 1.0e-100 # time increment between outputs + + +nlim = -1 # cycle limit +tlim = 100.0 # time limit (to be reset by problem generator) +integrator = rk1 # time integration algorithm +ncycle_out = 1 # interval for stdout summary info + + +# do_coalesced_comms = true +refinement = adaptive +numlevel = 3 +multigrid = true + +nx1 = 64 +x1min = -32.0 +x1max = 32.0 +ix1_bc = outflow +ox1_bc = outflow + +nx2 = 64 +x2min = -32.0 +x2max = 32.0 +ix2_bc = outflow +ox2_bc = outflow + +nx3 = 64 +x3min = -32.0 +x3max = 32.0 +ix3_bc = outflow +ox3_bc = outflow + + +nx1 = 16 +nx2 = 16 +nx3 = 16 + + +method = magnitude +comparator = greater_than +field = gas.prim.density_0 +refine_tol = 0.1 +derefine_tol = 0.1 + + +gas = true +gravity = true + + +cfl = 0.9 +reconstruct = plm +riemann = hllc +gamma = 2.0 + + +precondition = true +print_per_step = true +max_iterations = 100 +residual_tolerance = 1.0e-20 +smoother = SRJ2 +do_FAS = true +units_override = true +use_swindle = false +ix1_bc = zero +ox1_bc = zero +ix2_bc = zero +ox2_bc = zero +ix3_bc = zero +ox3_bc = zero + + +iprob = 1 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 959c9ff7..a7ed7dc0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,10 +40,13 @@ set (SRC_LIST gravity/binary_mass.cpp gravity/gravity.cpp + gravity/gravity.hpp gravity/nbody_gravity.hpp gravity/point_mass.cpp + gravity/poisson_driver.cpp + gravity/poisson_equation.hpp + gravity/self_gravity.cpp gravity/uniform.cpp - gravity/gravity.hpp nbody/nbody.cpp nbody/nbody.hpp @@ -64,6 +67,7 @@ set (SRC_LIST pgen/kh.hpp pgen/linear_wave.hpp pgen/lw.hpp + pgen/polytrope.hpp pgen/pgen.hpp pgen/problem_modifier.hpp pgen/rt.hpp diff --git a/src/artemis.cpp b/src/artemis.cpp index 652ab1db..589414fa 100644 --- a/src/artemis.cpp +++ b/src/artemis.cpp @@ -80,6 +80,7 @@ Packages_t ProcessPackages(std::unique_ptr &pin) { artemis->AddParam("x2max", x2max); artemis->AddParam("x3min", x3min); artemis->AddParam("x3max", x3max); + artemis->AddParam("domain_volume", (x1max - x1min) * (x2max - x2min) * (x3max - x3min)); // Add optionally enrollable operator split Metadata flag parthenon::MetadataFlag MetadataOperatorSplit = diff --git a/src/artemis.hpp b/src/artemis.hpp index 0fd8cf72..46942ab9 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -85,6 +85,11 @@ ARTEMIS_VARIABLE(rad.opac, scattering); } // namespace opac } // namespace rad +namespace grav { +ARTEMIS_VARIABLE(grav, phi); +ARTEMIS_VARIABLE(grav, rhs); +} // namespace grav + #undef ARTEMIS_VARIABLE // Restart options (see Parthenon #1231) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index 2cce1682..0b98623d 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -115,8 +115,10 @@ TaskListStatus ArtemisDriver::Step() { // Prepare registers PreStepTasks(); + TaskListStatus status; + // Execute explicit, unsplit physics - auto status = StepTasks().Execute(); + status = StepTasks().Execute(); if (status != TaskListStatus::complete) return status; // Operator split, background linear advection (for shearing box) @@ -206,6 +208,9 @@ TaskCollection ArtemisDriver::StepTasks() { const Real g1 = integrator->gam1[stage - 1]; const Real bdt = integrator->beta[stage - 1] * integrator->dt; + // Compute gravitational potential + if (do_gravity) Gravity::SolvePoisson(tc, pmesh); + TaskRegion &tr = tc.AddRegion(num_partitions); for (int i = 0; i < num_partitions; i++) { auto &tl = tr[i]; diff --git a/src/gravity/gravity.cpp b/src/gravity/gravity.cpp index 03aaa9bd..d0e1cb14 100644 --- a/src/gravity/gravity.cpp +++ b/src/gravity/gravity.cpp @@ -11,14 +11,43 @@ // the public, perform publicly and display publicly, and to permit others to do so. //======================================================================================== +// Parthenon includes +#include +#include +#include +#include +#include +#include +#include +#include + // Artemis includes -#include "gravity/gravity.hpp" #include "artemis.hpp" #include "geometry/geometry.hpp" +#include "gravity/gravity.hpp" #include "gravity/nbody_gravity.hpp" +#include "gravity/poisson_equation.hpp" namespace Gravity { +using namespace parthenon::BoundaryFunction; + +struct any_poisson : public parthenon::variable_names::base_t { + template + KOKKOS_INLINE_FUNCTION any_poisson(Ts &&...args) + : base_t(std::forward(args)...) {} + static std::string name() { return "grav[.].*"; } +}; + +template +auto Zero() { + return [](std::shared_ptr> &rc, bool coarse) -> void { + using namespace parthenon; + using namespace parthenon::BoundaryFunction; + GenericBC(rc, coarse, 0.0); + }; +} + //---------------------------------------------------------------------------------------- //! \fn StateDescriptor Gravity::Initialize //! \brief Adds intialization function for gravity package @@ -118,9 +147,89 @@ std::shared_ptr Initialize(ParameterInput *pin, auto &nbody_pkg = packages.Get("nbody"); params.Add("gm", nbody_pkg->Param("gm")); } + if (pin->DoesBlockExist("gravity/self")) { + count++; + gtype = GravityType::self; + block_name = "gravity/self"; + PARTHENON_REQUIRE(coords == Coordinates::cartesian, + "Self-gravity currently only supports Cartesian coordinates"); - PARTHENON_REQUIRE((count > 0) && (gtype != GravityType::null), "Unknown gravity node!"); + // Units and constants + const Real gcode = constants.GetGCode(); + Real four_pi_G = 4.0 * M_PI * gcode; + if (pin->GetOrAddBoolean(block_name, "units_override", false)) { + four_pi_G = pin->GetOrAddReal(block_name, "four_pi_G", 1.0); + } + params.Add("four_pi_G", four_pi_G); + + // Jeans swindle + bool pi1 = (pin->GetOrAddString("parthenon/mesh", "ix1_bc", "outflow") == "periodic"); + bool po1 = (pin->GetOrAddString("parthenon/mesh", "ox1_bc", "outflow") == "periodic"); + bool pi2 = (pin->GetOrAddString("parthenon/mesh", "ix2_bc", "outflow") == "periodic"); + bool po2 = (pin->GetOrAddString("parthenon/mesh", "ox2_bc", "outflow") == "periodic"); + bool pi3 = (pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow") == "periodic"); + bool po3 = (pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow") == "periodic"); + const bool needs_swindle = (pi1 && pi2 && pi3 && po1 && po2 && po3); + const bool swindle = pin->GetOrAddBoolean(block_name, "use_swindle", needs_swindle); + params.Add("use_swindle", swindle); + PARTHENON_REQUIRE(swindle || !(needs_swindle), + "Fully periodic BCs require Jeans swindle!"); + if (swindle && !(needs_swindle)) { + PARTHENON_WARN("Invoking Jeans swindle when not mandated by BCs"); + } + + // Gravity Package FillDerived function + gravity->FillDerivedMesh = FillPoissonRHS; + + // Enroll ZeroBC + using BF = parthenon::BoundaryFace; + constexpr auto IN = BCSide::Inner; + constexpr auto OUT = BCSide::Outer; + const bool zi1 = (pin->GetOrAddString(block_name, "ix1_bc", "default") == "zero"); + const bool zo1 = (pin->GetOrAddString(block_name, "ox1_bc", "default") == "zero"); + const bool zi2 = (pin->GetOrAddString(block_name, "ix2_bc", "default") == "zero"); + const bool zo2 = (pin->GetOrAddString(block_name, "ox2_bc", "default") == "zero"); + const bool zi3 = (pin->GetOrAddString(block_name, "ix3_bc", "default") == "zero"); + const bool zo3 = (pin->GetOrAddString(block_name, "ox3_bc", "default") == "zero"); + if (zi1) gravity->UserBoundaryFunctions[BF::inner_x1].push_back(Zero()); + if (zo1) gravity->UserBoundaryFunctions[BF::inner_x2].push_back(Zero()); + if (zi2) gravity->UserBoundaryFunctions[BF::inner_x3].push_back(Zero()); + if (zo2) gravity->UserBoundaryFunctions[BF::outer_x1].push_back(Zero()); + if (zi3) gravity->UserBoundaryFunctions[BF::outer_x2].push_back(Zero()); + if (zo3) gravity->UserBoundaryFunctions[BF::outer_x3].push_back(Zero()); + + // Gravitational potential + using namespace parthenon::refinement_ops; + std::vector flags{Metadata::Cell, Metadata::Independent, + Metadata::FillGhost, Metadata::WithFluxes, + Metadata::GMGRestrict, Metadata::GMGProlongate}; + Metadata m = Metadata(flags); + m.RegisterRefinementOps(); + gravity->AddField(m); + + // 4piG * \Sum rho + auto mrhs = Metadata( + {Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::FillGhost}); + gravity->AddField(mrhs); + + // Multigrid + using PoissEq = PoissonEquation; + PoissEq eq(pin, "poisson"); + params.Add("poisson_equation", eq, parthenon::Params::Mutability::Mutable); + std::shared_ptr psolver; + using prolongator_t = parthenon::solvers::ProlongationBlockInteriorZeroDirichlet; + using preconditioner_t = parthenon::solvers::MGSolver; + psolver = + std::make_shared>( + "base", "phi", "rhs", pin, block_name, PoissEq(pin, block_name)); + // psolver = std::make_shared>( + // "base", "phi", "rhs", pin, block_name, PoissEq(pin, block_name)); + + params.Add("solver_pointer", psolver); + } + + PARTHENON_REQUIRE((count > 0) && (gtype != GravityType::null), "Unknown gravity node!"); PARTHENON_REQUIRE(count == 1, "artemis only supports 1 gravity type at this time"); params.Add("type", gtype); @@ -155,11 +264,91 @@ TaskStatus ExternalGravity(MeshData *md, const Real time, const Real dt) { if (pkg->template Param("npart") > 0) return NBodyGravity(md, time, dt); } + } else if (gtype == GravityType::self) { + return SelfGravity(md, time, dt); } } return TaskStatus::complete; } +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus Gravity::ExternalGravity +//! \brief Wrapper function for external gravity options +template +void FillPoissonRHS(MeshData *md) { + using parthenon::MakePackDescriptor; + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + + // Extract artemis parameters + 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"); + + // Packing and indexing + static auto desc = + MakePackDescriptor( + resolved_pkgs.get()); + auto vmesh = desc.GetPack(md); + const int nblocks = md->NumBlocks(); + IndexRange ib = md->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBoundsK(IndexDomain::interior); + IndexRange ibe = md->GetBoundsI(IndexDomain::entire); + IndexRange jbe = md->GetBoundsJ(IndexDomain::entire); + IndexRange kbe = md->GetBoundsK(IndexDomain::entire); + + // Compute mean density + Real grav_mean_rho = 0.0; + auto &grav_pkg = pm->packages.Get("gravity"); + const bool use_swindle = grav_pkg->template Param("use_swindle"); + if (use_swindle) { + Real total_mass = 0.0; + Real total_volume = artemis_pkg->template Param("domain_volume"); + parthenon::par_reduce( + parthenon::loop_pattern_mdrange_tag, "Gravity::TotalMass", + parthenon::DevExecSpace(), 0, nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i, + Real &mtot) { + geometry::Coords coords(cpars, vmesh.GetCoordinates(b), k, j, i); + const Real vv = coords.Volume(); + for (int n = 0; n < do_gas * vmesh.GetSize(b, gas::prim::density()); ++n) { + mtot += vmesh(b, gas::prim::density(n), k, j, i) * vv; + } + for (int n = 0; n < do_dust * vmesh.GetSize(b, dust::prim::density()); ++n) { + mtot += vmesh(b, dust::prim::density(n), k, j, i) * vv; + } + }, + Kokkos::Sum(total_mass)); + Kokkos::fence(); + +#ifdef MPI_PARALLEL + // NOTE(@pdmullen): This reduction only works because we require pack_size==-1... + MPI_Allreduce(MPI_IN_PLACE, &total_mass, 1, MPI_PARTHENON_REAL, MPI_SUM, + MPI_COMM_WORLD); +#endif + grav_mean_rho = total_mass / total_volume; + } + + // Set RHS + const Real four_pi_G = grav_pkg->template Param("four_pi_G"); + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SetRHS", parthenon::DevExecSpace(), 0, nblocks - 1, kbe.s, + kbe.e, jbe.s, jbe.e, ibe.s, ibe.e, + KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) { + Real &rhs = vmesh(b, grav::rhs(), k, j, i) = 0.0; + for (int n = 0; n < do_gas * vmesh.GetSize(b, gas::prim::density()); ++n) { + rhs += vmesh(b, gas::prim::density(n), k, j, i); + } + for (int n = 0; n < do_dust * vmesh.GetSize(b, dust::prim::density()); ++n) { + rhs += vmesh(b, dust::prim::density(n), k, j, i); + } + rhs -= use_swindle * grav_mean_rho; + rhs *= four_pi_G; + }); +} + //---------------------------------------------------------------------------------------- //! template instantiations typedef Coordinates G; @@ -170,5 +359,11 @@ template TaskStatus ExternalGravity(MD *m, const Real t, const R template TaskStatus ExternalGravity(MD *m, const Real t, const Real dt); template TaskStatus ExternalGravity(MD *m, const Real t, const Real dt); template TaskStatus ExternalGravity(MD *m, const Real t, const Real dt); +template void FillPoissonRHS(MD *m); +template void FillPoissonRHS(MD *m); +template void FillPoissonRHS(MD *m); +template void FillPoissonRHS(MD *m); +template void FillPoissonRHS(MD *m); +template void FillPoissonRHS(MD *m); } // namespace Gravity diff --git a/src/gravity/gravity.hpp b/src/gravity/gravity.hpp index 872ad44e..ba112b84 100644 --- a/src/gravity/gravity.hpp +++ b/src/gravity/gravity.hpp @@ -22,7 +22,7 @@ namespace Gravity { -enum class GravityType { uniform, point, binary, nbody, null }; +enum class GravityType { uniform, point, binary, nbody, self, null }; //---------------------------------------------------------------------------------------- //! \struct Orbit @@ -113,8 +113,15 @@ TaskStatus NBodyGravityFixed(MeshData *md, const Real time, const Real dt) template TaskStatus ExternalGravity(MeshData *md, const Real time, const Real dt); -KOKKOS_INLINE_FUNCTION -Real quad_ramp(const Real x) { return SQR(x); } +template +void FillPoissonRHS(MeshData *md); + +template +TaskStatus SelfGravity(MeshData *md, const Real time, const Real dt); + +void SolvePoisson(TaskCollection &tc, Mesh *pmesh); + +KOKKOS_INLINE_FUNCTION Real quad_ramp(const Real x) { return SQR(x); } } // namespace Gravity diff --git a/src/gravity/poisson_driver.cpp b/src/gravity/poisson_driver.cpp new file mode 100644 index 00000000..257204a0 --- /dev/null +++ b/src/gravity/poisson_driver.cpp @@ -0,0 +1,77 @@ +//======================================================================================== +// (C) (or copyright) 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 +#include +#include +#include +#include + +// Parthenon includes +#include +#include +#include +#include +#include +#include +#include +#include + +// Artemis includes +#include "gravity/gravity.hpp" +#include "gravity/poisson_equation.hpp" + +using namespace parthenon::driver::prelude; + +namespace Gravity { + +//---------------------------------------------------------------------------------------- +//! \fn TaskListStatus Gravity::PoissonDriver +//! \brief +void SolvePoisson(TaskCollection &tc, Mesh *pmesh) { + using namespace parthenon; + TaskID none(0); + + auto pkg = pmesh->packages.Get("gravity"); + auto psolver = + pkg->Param>("solver_pointer"); + + auto partitions = pmesh->GetDefaultBlockPartitions(); + const int num_partitions = partitions.size(); + TaskRegion ®ion = tc.AddRegion(num_partitions); + for (int i = 0; i < num_partitions; ++i) { + TaskList &tl = region[i]; + auto &md = pmesh->mesh_data.Add("base", partitions[i]); + auto &md_phi = pmesh->mesh_data.Add("phi", md, {grav::phi::name()}); + auto &md_rhs = pmesh->mesh_data.Add("rhs", md, {grav::phi::name()}); + + // Move the rhs variable into the rhs stage for stage based solver + auto copy_rhs = tl.AddTask( + none, TF(solvers::utils::between_fields::CopyData), md); + copy_rhs = + tl.AddTask(copy_rhs, TF(solvers::utils::CopyData>), + md, md_rhs); + + // Set initial solution guess to zero + auto zero_phi = + tl.AddTask(copy_rhs, TF(solvers::utils::SetToZero), md_phi); + auto setup = psolver->AddSetupTasks(tl, zero_phi, i, pmesh); + auto solve = psolver->AddTasks(tl, setup, i, pmesh); + + // Move the solution back so it is output + auto copy_back = tl.AddTask( + solve, TF(solvers::utils::CopyData>), md_phi, md); + } +} + +} // namespace Gravity diff --git a/src/gravity/poisson_equation.hpp b/src/gravity/poisson_equation.hpp new file mode 100644 index 00000000..6e6905b3 --- /dev/null +++ b/src/gravity/poisson_equation.hpp @@ -0,0 +1,222 @@ +//======================================================================================== +// (C) (or copyright) 2023-2024. 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 GRAVITY_POISSON_EQUATION_HPP_ +#define GRAVITY_POISSON_EQUATION_HPP_ + +// Parthenon includes +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Gravity { +constexpr parthenon::TopologicalElement te = parthenon::TopologicalElement::CC; + +// This class implement methods for calculating A.x = y and returning the diagonal of A, +// where A is the the matrix representing the discretized Poisson equation on the grid. +template +class PoissonEquation { + public: + using IndependentVars = parthenon::TypeList; + + PoissonEquation(parthenon::ParameterInput *pin, const std::string &label) {} + + // Add tasks to calculate the result of the matrix A (which is implicitly defined by + // this class) being applied to x_t and store it in field out_t + parthenon::TaskID Ax(parthenon::TaskList &tl, parthenon::TaskID depends_on, + std::shared_ptr> &md_mat, + std::shared_ptr> &md_in, + std::shared_ptr> &md_out) { + auto flux_res = tl.AddTask(depends_on, CalculateFluxes, md_mat, md_in); + if (!(md_mat->grid.type == parthenon::GridType::two_level_composite)) { + auto start_flxcor = + tl.AddTask(flux_res, parthenon::StartReceiveFluxCorrections, md_in); + auto send_flxcor = + tl.AddTask(flux_res, parthenon::LoadAndSendFluxCorrections, md_in); + auto recv_flxcor = + tl.AddTask(start_flxcor, parthenon::ReceiveFluxCorrections, md_in); + flux_res = tl.AddTask(recv_flxcor, parthenon::SetFluxCorrections, md_in); + } + return tl.AddTask(flux_res, FluxMultiplyMatrix, md_in, md_out); + } + + template + KOKKOS_INLINE_FUNCTION auto GetEffectiveInverseDx2(const coords_t &coords, const int k, + const int j, const int i) { + using TE = parthenon::TopologicalElement; + constexpr TE te = dir == X1DIR ? TE::F1 : (dir == X2DIR ? TE::F2 : TE::F3); + constexpr int ioff = (dir == X1DIR); + constexpr int joff = (dir == X2DIR); + constexpr int koff = (dir == X3DIR); + + const Real xp = coords.template Xc(k + koff, j + joff, i + ioff); + const Real xc = coords.template Xc(k, j, i); + const Real xm = coords.template Xc(k - koff, j - joff, i - ioff); + + const Real dxp = xp - xc; + const Real dxm = xc - xm; + const Real Ap = coords.template Volume(k + koff, j + joff, i + ioff); + const Real Am = coords.template Volume(k, j, i); + const Real Vol = coords.template Volume(k, j, i); + return std::make_pair(Ap / (dxp * Vol), Am / (dxm * Vol)); + } + + // Calculate an approximation to the diagonal of the matrix A and store it in diag_t. + // For a uniform grid or when flux correction is ignored, this diagonal calculation + // is exact. Exactness is (probably) not required since it is just used in Jacobi + // iterations. + parthenon::TaskStatus SetDiagonal(std::shared_ptr> &md_mat, + std::shared_ptr> &md_diag) { + using namespace parthenon; + const int ndim = md_mat->GetMeshPointer()->ndim; + IndexRange ib = md_mat->GetBoundsI(IndexDomain::interior, te); + IndexRange jb = md_mat->GetBoundsJ(IndexDomain::interior, te); + IndexRange kb = md_mat->GetBoundsK(IndexDomain::interior, te); + int nblocks = md_mat->NumBlocks(); + + auto desc_diag = parthenon::MakePackDescriptor(md_diag.get()); + auto pack_diag = desc_diag.GetPack(md_diag.get()); + using TE = parthenon::TopologicalElement; + parthenon::par_for( + "StoreDiagonal", 0, pack_diag.GetNBlocks() - 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) { + const auto &coords = pack_diag.GetCoordinates(b); + // Build the unigrid diagonal of the matrix + Real diag_elem = 0.0; + { + auto [idx2p, idx2m] = GetEffectiveInverseDx2(coords, k, j, i); + diag_elem -= (idx2m + idx2p); + } + if (ndim > 1) { + auto [idx2p, idx2m] = GetEffectiveInverseDx2(coords, k, j, i); + diag_elem -= (idx2m + idx2p); + } + if (ndim > 2) { + auto [idx2p, idx2m] = GetEffectiveInverseDx2(coords, k, j, i); + diag_elem -= (idx2m + idx2p); + } + pack_diag(b, te, var_t(), k, j, i) = diag_elem; + }); + return TaskStatus::complete; + } + + static parthenon::TaskStatus + CalculateFluxes(std::shared_ptr> &md_mat, + std::shared_ptr> &md) { + using namespace parthenon; + const int ndim = md->GetMeshPointer()->ndim; + using TE = parthenon::TopologicalElement; + TE te = TE::CC; + IndexRange ib = md->GetBoundsI(IndexDomain::interior, te); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior, te); + IndexRange kb = md->GetBoundsK(IndexDomain::interior, te); + int nblocks = md->NumBlocks(); + + auto desc = parthenon::MakePackDescriptor(md.get(), {}, {PDOpt::WithFluxes}); + auto pack = desc.GetPack(md.get()); + parthenon::par_for( + "CaclulateFluxes", 0, pack.GetNBlocks() - 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) { + const auto &coords = pack.GetCoordinates(b); + pack.flux(b, X1DIR, var_t(), k, j, i) = + (pack(b, te, var_t(), k, j, i - 1) - pack(b, te, var_t(), k, j, i)) / + coords.template Dxc(k, j, i); + if (i == ib.e) + pack.flux(b, X1DIR, var_t(), k, j, i + 1) = + (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k, j, i + 1)) / + coords.template Dxc(k, j, i + 1); + + if (ndim > 1) { + pack.flux(b, X2DIR, var_t(), k, j, i) = + (pack(b, te, var_t(), k, j - 1, i) - pack(b, te, var_t(), k, j, i)) / + coords.template Dxc(k, j, i); + if (j == jb.e) + pack.flux(b, X2DIR, var_t(), k, j + 1, i) = + (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k, j + 1, i)) / + coords.template Dxc(k, j + 1, i); + } + + if (ndim > 2) { + pack.flux(b, X3DIR, var_t(), k, j, i) = + (pack(b, te, var_t(), k - 1, j, i) - pack(b, te, var_t(), k, j, i)) / + coords.template Dxc(k, j, i); + if (k == kb.e) + pack.flux(b, X3DIR, var_t(), k + 1, j, i) = + (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k + 1, j, i)) / + coords.template Dxc(k + 1, j, i); + } + }); + return TaskStatus::complete; + } + + // Calculate A in_t = out_t (in the region covered by md) for a given set of fluxes + // calculated with in_t (which have possibly been corrected at coarse fine boundaries) + static parthenon::TaskStatus + FluxMultiplyMatrix(std::shared_ptr> &md, + std::shared_ptr> &md_out) { + using namespace parthenon; + const int ndim = md->GetMeshPointer()->ndim; + using TE = parthenon::TopologicalElement; + TE te = TE::CC; + IndexRange ib = md->GetBoundsI(IndexDomain::interior, te); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior, te); + IndexRange kb = md->GetBoundsK(IndexDomain::interior, te); + int nblocks = md->NumBlocks(); + + static auto desc = + parthenon::MakePackDescriptor(md.get(), {}, {PDOpt::WithFluxes}); + static auto desc_out = parthenon::MakePackDescriptor(md_out.get()); + auto pack = desc.GetPack(md.get()); + auto pack_out = desc_out.GetPack(md_out.get()); + parthenon::par_for( + "FluxMultiplyMatrix", 0, pack.GetNBlocks() - 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) { + const auto &coords = pack.GetCoordinates(b); + Real dx1 = coords.template Dxc(k, j, i); + pack_out(b, te, var_t(), k, j, i) = 0.0; + pack_out(b, te, var_t(), k, j, i) += + (pack.flux(b, X1DIR, var_t(), k, j, i) * + coords.template Volume(k, j, i) - + pack.flux(b, X1DIR, var_t(), k, j, i + 1) * + coords.template Volume(k, j, i + 1)) / + coords.template Volume(k, j, i); + + if (ndim > 1) { + pack_out(b, te, var_t(), k, j, i) += + (pack.flux(b, X2DIR, var_t(), k, j, i) * + coords.template Volume(k, j, i) - + pack.flux(b, X2DIR, var_t(), k, j + 1, i) * + coords.template Volume(k, j + 1, i)) / + coords.template Volume(k, j, i); + } + + if (ndim > 2) { + pack_out(b, te, var_t(), k, j, i) += + (pack.flux(b, X3DIR, var_t(), k, j, i) * + coords.template Volume(k, j, i) - + pack.flux(b, X3DIR, var_t(), k + 1, j, i) * + coords.template Volume(k + 1, j, i)) / + coords.template Volume(k, j, i); + } + }); + return TaskStatus::complete; + } +}; + +} // namespace Gravity + +#endif // GRAVITY_POISSON_EQUATION_HPP_ diff --git a/src/gravity/self_gravity.cpp b/src/gravity/self_gravity.cpp new file mode 100644 index 00000000..2426a070 --- /dev/null +++ b/src/gravity/self_gravity.cpp @@ -0,0 +1,145 @@ +//======================================================================================== +// (C) (or copyright) 2023-2025. 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. +//======================================================================================== + +// Artemis includes +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "gravity/gravity.hpp" +#include "utils/artemis_utils.hpp" + +using namespace parthenon::package::prelude; +using ArtemisUtils::VI; + +namespace Gravity { +//---------------------------------------------------------------------------------------- +//! \fn TaskStatus Gravity::SelfGravity +//! \brief Applies accelerations due to a constant g +template +TaskStatus SelfGravity(MeshData *md, const Real time, const Real dt) { + 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()); + static auto descf = MakePackDescriptor( + resolved_pkgs.get(), {}, {parthenon::PDOpt::WithFluxes}); + auto vmesh = desc.GetPack(md); + auto vflux = descf.GetPack(md); + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + const bool multi_d = (pm->ndim > 1); + const bool three_d = (pm->ndim > 2); + const int d1 = X1DIR; + const int d2 = d1 + multi_d; + const int d3 = d2 + three_d; + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "SelfGravity", 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 &dx = coords.GetCellWidths(); + const Real hdtodx1 = (0.5 * dt / dx[0]); + const Real hdtodx2 = multi_d * (0.5 * dt / dx[1]); + const Real hdtodx3 = three_d * (0.5 * dt / dx[2]); + + // Potential differences + const Real &phic = vmesh(b, grav::phi(), k, j, i); + const Real dpl1 = -(phic - vmesh(b, grav::phi(), k, j, i - 1)); + const Real dpr1 = -(vmesh(b, grav::phi(), k, j, i + 1) - phic); + const Real dpl2 = -(phic - vmesh(b, grav::phi(), k, j - multi_d, i)); + const Real dpr2 = -(vmesh(b, grav::phi(), k, j + multi_d, i) - phic); + const Real dpl3 = -(phic - vmesh(b, grav::phi(), k - three_d, j, i)); + const Real dpr3 = -(vmesh(b, grav::phi(), k + three_d, j, i) - phic); + + // Shared b/w gas and dust + const Real wdt1 = hdtodx1 * (dpl1 + dpr1); + const Real wdt2 = hdtodx2 * (dpl2 + dpr2); + const Real wdt3 = hdtodx3 * (dpl3 + dpr3); + + // Check answer + const Real rho_alt = + (6.0 * vmesh(b, grav::phi(), k, j, i) - vmesh(b, grav::phi(), k, j, i - 1) - + vmesh(b, grav::phi(), k, j, i + 1) - + vmesh(b, grav::phi(), k, j - multi_d, i) - + vmesh(b, grav::phi(), k, j + multi_d, i) - + vmesh(b, grav::phi(), k - three_d, j, i) - + vmesh(b, grav::phi(), k + three_d, j, i)) / + (SQR(dx[0])); + const Real my_rho = vmesh(b, gas::prim::density(0), k, j, i); + if (std::abs(my_rho - rho_alt) / rho_alt > 1.0e-5) { + printf("%d %d %d %24.16e %24.16e %24.16e\n", k, j, i, rho_alt, my_rho, + vmesh(b, grav::phi(), k, j, i)); + std::exit(1); + } + + // if (do_gas) { + // // Gravitational acceleration and energy release + // for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { + // const Real &rr = vmesh(b, gas::prim::density(n), k, j, i); + // vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; + // vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; + // vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; + // vmesh(b, gas::cons::total_energy(n), k, j, i) += + // hdtodx1 * (vflux.flux(b, d1, gas::cons::density(n), k, j, i) * dpl1 + + // vflux.flux(b, d1, gas::cons::density(n), k, j, i + 1) * + // dpr1); + // vmesh(b, gas::cons::total_energy(n), k, j, i) += + // hdtodx2 * + // (vflux.flux(b, d2, gas::cons::density(n), k, j, i) * dpl2 + + // vflux.flux(b, d2, gas::cons::density(n), k, j + multi_d, i) * dpr2); + // vmesh(b, gas::cons::total_energy(n), k, j, i) += + // hdtodx3 * + // (vflux.flux(b, d3, gas::cons::density(n), k, j, i) * dpl3 + + // vflux.flux(b, d3, gas::cons::density(n), k + three_d, j, i) * dpr3); + // } + // } + + // if (do_dust) { + // for (int n = 0; n < vmesh.GetSize(b, dust::prim::density()); ++n) { + // // Gravitational acceleration + // const Real &rr = vmesh(b, dust::prim::density(n), k, j, i); + // vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; + // vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; + // vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; + // } + // } + }); + + printf("Success!\n"); + + return TaskStatus::complete; +} + +//---------------------------------------------------------------------------------------- +//! template instantiations +typedef Coordinates G; +typedef MeshData MD; +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); +template TaskStatus SelfGravity(MD *m, const Real t, const Real d); + +} // namespace Gravity diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index d7240766..f82d045f 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -28,6 +28,7 @@ #include "kh.hpp" #include "linear_wave.hpp" #include "lw.hpp" +#include "polytrope.hpp" #include "rt.hpp" #include "shock.hpp" #include "strat.hpp" @@ -60,6 +61,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { linear_wave::ProblemGenerator(pmb, pin); } else if (name == "lw") { lw::ProblemGenerator(pmb, pin); + } else if (name == "polytrope") { + polytrope::ProblemGenerator(pmb, pin); } else if (name == "kh") { kh::ProblemGenerator(pmb, pin); } else if (name == "rt") { diff --git a/src/pgen/polytrope.hpp b/src/pgen/polytrope.hpp new file mode 100644 index 00000000..c3ab2075 --- /dev/null +++ b/src/pgen/polytrope.hpp @@ -0,0 +1,113 @@ +//======================================================================================== +// (C) (or copyright) 2023-2025. 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 PGEN_POLYTROPE_HPP_ +#define PGEN_POLYTROPE_HPP_ +//! \file polytrope.hpp + +// C/C++ headers +#include +#include +#include +#include +#include +#include +#include + +// Artemis headers +#include "artemis.hpp" +#include "geometry/geometry.hpp" +#include "utils/artemis_utils.hpp" + +namespace polytrope { + +//---------------------------------------------------------------------------------------- +//! \fn void ProblemGenerator::LinearWave_() +//! \brief Sets initial conditions for polytrope tests +template +inline void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { + using parthenon::MakePackDescriptor; + const Mesh *pmesh = pmb->pmy_mesh; + const int ndim = pmesh->ndim; + using TE = parthenon::TopologicalElement; + + // packing and capture variables for kernel + auto &md = pmb->meshblock_data.Get(); + for (auto &var : md->GetVariableVector()) { + if (!var->IsAllocated()) pmb->AllocateSparse(var->label()); + } + static auto desc = + MakePackDescriptor( + (pmb->resolved_packages).get()); + auto v = desc.GetPack(md.get()); + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + auto &pco = pmb->coords; + const auto &cpars = + pmb->packages.Get("artemis")->template Param("coord_params"); + + // Polytrope params + const int iprob = pin->GetOrAddInteger("problem", "iprob", 1); + const Real x10a = pin->GetOrAddReal("problem", "x10a", (iprob == 2) * 4.0); + const Real x20a = pin->GetOrAddReal("problem", "x20a", (iprob == 2) * 2.5); + const Real x30a = pin->GetOrAddReal("problem", "x30a", (iprob == 2) * 0.0); + const Real x10b = pin->GetOrAddReal("problem", "x10b", -4.0); + const Real x20b = pin->GetOrAddReal("problem", "x20b", -2.5); + const Real x30b = pin->GetOrAddReal("problem", "x30b", 0.0); + const Real rho_amb = pin->GetOrAddReal("problem", "rho_amb", 1.0e-3); + const Real sie_amb = pin->GetOrAddReal("problem", "sie_amb", 1.0e2); + + // NOTE(@pdmullen): Hardcoded params related to n=1 polytrope profile... + const Real alpha = std::sqrt(0.5); + const Real cutoff = 0.75 * M_PI; + + // Polytrope init + pmb->par_for( + "polytrope", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int k, const int j, const int i) { + // cell-centered coordinates + geometry::Coords coords(cpars, pco, k, j, i); + const auto &xv = coords.GetCellCenter(); + + // Compute Lane-Emden profiles + const Real ar1 = + alpha * std::sqrt(SQR(xv[0] - x10a) + SQR(xv[1] - x20a) + SQR(xv[2] - x30a)); + const Real ar2 = + alpha * std::sqrt(SQR(xv[0] - x10b) + SQR(xv[1] - x20b) + SQR(xv[2] - x30b)); + const Real lane_emden1 = std::sin(ar1) / ar1; + const Real lane_emden2 = std::sin(ar2) / ar2; + const bool inside1 = (ar1 < cutoff); + const bool inside2 = (ar2 < cutoff); + + // Initialize polytrope(s) + Real trho = Null(), tsie = Null(); + if (iprob == 1) { + trho = inside1 ? lane_emden1 : rho_amb; + tsie = inside1 ? lane_emden1 : sie_amb; + } else if (iprob == 2) { + trho = inside1 ? lane_emden1 : (inside2 ? lane_emden2 : rho_amb); + tsie = inside1 ? lane_emden1 : (inside2 ? lane_emden2 : sie_amb); + } else { + PARTHENON_FAIL("iprob not recognized!"); + } + v(0, gas::prim::density(), k, j, i) = trho; + v(0, gas::prim::sie(), k, j, i) = tsie; + v(0, gas::prim::velocity(0), k, j, i) = (iprob == 2) * (inside2 - inside1); + v(0, gas::prim::velocity(1), k, j, i) = 0.0; + v(0, gas::prim::velocity(2), k, j, i) = 0.0; + }); +} + +} // namespace polytrope + +#endif // PGEN_POLYTROPE_HPP_ diff --git a/src/utils/integrators/artemis_integrator.hpp b/src/utils/integrators/artemis_integrator.hpp index 290cc942..2607e2ce 100644 --- a/src/utils/integrators/artemis_integrator.hpp +++ b/src/utils/integrators/artemis_integrator.hpp @@ -88,19 +88,19 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const Real g0, const Real bdt_vol = beta_dt / coords.Volume(); // Advance state vector with flux divergence - for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { - Real &v0n = v0(b, n, k, j, i); - Real &v1n = v1(b, n, k, j, i); - v0n = g0 * v0n + g1 * v1n; - if constexpr (include_divf) { - v0n += bdt_vol * ((ax1[0] * v0.flux(b, d1, n, k, j, i) - - ax1[1] * v0.flux(b, d1, n, k, j, i + 1)) + - (ax2[0] * v0.flux(b, d2, n, k, j, i) - - ax2[1] * v0.flux(b, d2, n, k, j + multi_d, i)) + - (ax3[0] * v0.flux(b, d3, n, k, j, i) - - ax3[1] * v0.flux(b, d3, n, k + three_d, j, i))); - } - } + // for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { + // Real &v0n = v0(b, n, k, j, i); + // Real &v1n = v1(b, n, k, j, i); + // v0n = g0 * v0n + g1 * v1n; + // if constexpr (include_divf) { + // v0n += bdt_vol * ((ax1[0] * v0.flux(b, d1, n, k, j, i) - + // ax1[1] * v0.flux(b, d1, n, k, j, i + 1)) + + // (ax2[0] * v0.flux(b, d2, n, k, j, i) - + // ax2[1] * v0.flux(b, d2, n, k, j + multi_d, i)) + + // (ax3[0] * v0.flux(b, d3, n, k, j, i) - + // ax3[1] * v0.flux(b, d3, n, k + three_d, j, i))); + // } + // } }); return TaskStatus::complete; } From e2659cee163829456fc3e95b12d9e84b378043f6 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 17 Dec 2025 13:51:31 -0500 Subject: [PATCH 2/8] Add boundary exchange task after Poisson solve --- src/gravity/poisson_driver.cpp | 5 ++- src/gravity/self_gravity.cpp | 75 +++++++++++++--------------------- 2 files changed, 32 insertions(+), 48 deletions(-) diff --git a/src/gravity/poisson_driver.cpp b/src/gravity/poisson_driver.cpp index 257204a0..7be30296 100644 --- a/src/gravity/poisson_driver.cpp +++ b/src/gravity/poisson_driver.cpp @@ -68,9 +68,12 @@ void SolvePoisson(TaskCollection &tc, Mesh *pmesh) { auto setup = psolver->AddSetupTasks(tl, zero_phi, i, pmesh); auto solve = psolver->AddTasks(tl, setup, i, pmesh); + // Set BCs after solve + auto bcs = parthenon::AddBoundaryExchangeTasks(solve, tl, md_phi, pmesh->multilevel); + // Move the solution back so it is output auto copy_back = tl.AddTask( - solve, TF(solvers::utils::CopyData>), md_phi, md); + bcs, TF(solvers::utils::CopyData>), md_phi, md); } } diff --git a/src/gravity/self_gravity.cpp b/src/gravity/self_gravity.cpp index 2426a070..080a0da8 100644 --- a/src/gravity/self_gravity.cpp +++ b/src/gravity/self_gravity.cpp @@ -77,57 +77,38 @@ TaskStatus SelfGravity(MeshData *md, const Real time, const Real dt) { const Real wdt2 = hdtodx2 * (dpl2 + dpr2); const Real wdt3 = hdtodx3 * (dpl3 + dpr3); - // Check answer - const Real rho_alt = - (6.0 * vmesh(b, grav::phi(), k, j, i) - vmesh(b, grav::phi(), k, j, i - 1) - - vmesh(b, grav::phi(), k, j, i + 1) - - vmesh(b, grav::phi(), k, j - multi_d, i) - - vmesh(b, grav::phi(), k, j + multi_d, i) - - vmesh(b, grav::phi(), k - three_d, j, i) - - vmesh(b, grav::phi(), k + three_d, j, i)) / - (SQR(dx[0])); - const Real my_rho = vmesh(b, gas::prim::density(0), k, j, i); - if (std::abs(my_rho - rho_alt) / rho_alt > 1.0e-5) { - printf("%d %d %d %24.16e %24.16e %24.16e\n", k, j, i, rho_alt, my_rho, - vmesh(b, grav::phi(), k, j, i)); - std::exit(1); + if (do_gas) { + // Gravitational acceleration and energy release + for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { + const Real &rr = vmesh(b, gas::prim::density(n), k, j, i); + vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; + vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; + vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; + vmesh(b, gas::cons::total_energy(n), k, j, i) += + hdtodx1 * (vflux.flux(b, d1, gas::cons::density(n), k, j, i) * dpl1 + + vflux.flux(b, d1, gas::cons::density(n), k, j, i + 1) * dpr1); + vmesh(b, gas::cons::total_energy(n), k, j, i) += + hdtodx2 * + (vflux.flux(b, d2, gas::cons::density(n), k, j, i) * dpl2 + + vflux.flux(b, d2, gas::cons::density(n), k, j + multi_d, i) * dpr2); + vmesh(b, gas::cons::total_energy(n), k, j, i) += + hdtodx3 * + (vflux.flux(b, d3, gas::cons::density(n), k, j, i) * dpl3 + + vflux.flux(b, d3, gas::cons::density(n), k + three_d, j, i) * dpr3); + } } - // if (do_gas) { - // // Gravitational acceleration and energy release - // for (int n = 0; n < vmesh.GetSize(b, gas::prim::density()); ++n) { - // const Real &rr = vmesh(b, gas::prim::density(n), k, j, i); - // vmesh(b, gas::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; - // vmesh(b, gas::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; - // vmesh(b, gas::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; - // vmesh(b, gas::cons::total_energy(n), k, j, i) += - // hdtodx1 * (vflux.flux(b, d1, gas::cons::density(n), k, j, i) * dpl1 + - // vflux.flux(b, d1, gas::cons::density(n), k, j, i + 1) * - // dpr1); - // vmesh(b, gas::cons::total_energy(n), k, j, i) += - // hdtodx2 * - // (vflux.flux(b, d2, gas::cons::density(n), k, j, i) * dpl2 + - // vflux.flux(b, d2, gas::cons::density(n), k, j + multi_d, i) * dpr2); - // vmesh(b, gas::cons::total_energy(n), k, j, i) += - // hdtodx3 * - // (vflux.flux(b, d3, gas::cons::density(n), k, j, i) * dpl3 + - // vflux.flux(b, d3, gas::cons::density(n), k + three_d, j, i) * dpr3); - // } - // } - - // if (do_dust) { - // for (int n = 0; n < vmesh.GetSize(b, dust::prim::density()); ++n) { - // // Gravitational acceleration - // const Real &rr = vmesh(b, dust::prim::density(n), k, j, i); - // vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; - // vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; - // vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; - // } - // } + if (do_dust) { + for (int n = 0; n < vmesh.GetSize(b, dust::prim::density()); ++n) { + // Gravitational acceleration + const Real &rr = vmesh(b, dust::prim::density(n), k, j, i); + vmesh(b, dust::cons::momentum(VI(n, 0)), k, j, i) += rr * wdt1; + vmesh(b, dust::cons::momentum(VI(n, 1)), k, j, i) += rr * wdt2; + vmesh(b, dust::cons::momentum(VI(n, 2)), k, j, i) += rr * wdt3; + } + } }); - printf("Success!\n"); - return TaskStatus::complete; } From dd510a9880cb0d279556b55fd11311000001d0f9 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 17 Dec 2025 13:58:35 -0500 Subject: [PATCH 3/8] Cleanup --- inputs/self-gravity/poly.in | 4 +-- src/gravity/poisson_driver.cpp | 2 +- src/gravity/poisson_equation.hpp | 2 +- src/pgen/pgen.hpp | 2 +- src/utils/integrators/artemis_integrator.hpp | 26 ++++++++++---------- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/inputs/self-gravity/poly.in b/inputs/self-gravity/poly.in index 75ebf774..b24133b8 100644 --- a/inputs/self-gravity/poly.in +++ b/inputs/self-gravity/poly.in @@ -1,5 +1,5 @@ # ======================================================================================== -# (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved. +# (C) (or copyright) 2025. 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 @@ -12,7 +12,7 @@ # ======================================================================================== -problem = polytrope # name of the pgen +problem = polytrope # name of the pgen coordinates = cartesian # coordinate system diff --git a/src/gravity/poisson_driver.cpp b/src/gravity/poisson_driver.cpp index 7be30296..f4a5a480 100644 --- a/src/gravity/poisson_driver.cpp +++ b/src/gravity/poisson_driver.cpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2025. 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 diff --git a/src/gravity/poisson_equation.hpp b/src/gravity/poisson_equation.hpp index 6e6905b3..bf4f9a26 100644 --- a/src/gravity/poisson_equation.hpp +++ b/src/gravity/poisson_equation.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2025. 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 diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index f82d045f..35db513e 100644 --- a/src/pgen/pgen.hpp +++ b/src/pgen/pgen.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2025. 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 diff --git a/src/utils/integrators/artemis_integrator.hpp b/src/utils/integrators/artemis_integrator.hpp index 2607e2ce..290cc942 100644 --- a/src/utils/integrators/artemis_integrator.hpp +++ b/src/utils/integrators/artemis_integrator.hpp @@ -88,19 +88,19 @@ TaskStatus ApplyUpdate(MeshData *u0, MeshData *u1, const Real g0, const Real bdt_vol = beta_dt / coords.Volume(); // Advance state vector with flux divergence - // for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { - // Real &v0n = v0(b, n, k, j, i); - // Real &v1n = v1(b, n, k, j, i); - // v0n = g0 * v0n + g1 * v1n; - // if constexpr (include_divf) { - // v0n += bdt_vol * ((ax1[0] * v0.flux(b, d1, n, k, j, i) - - // ax1[1] * v0.flux(b, d1, n, k, j, i + 1)) + - // (ax2[0] * v0.flux(b, d2, n, k, j, i) - - // ax2[1] * v0.flux(b, d2, n, k, j + multi_d, i)) + - // (ax3[0] * v0.flux(b, d3, n, k, j, i) - - // ax3[1] * v0.flux(b, d3, n, k + three_d, j, i))); - // } - // } + for (int n = v0.GetLowerBound(b); n <= v0.GetUpperBound(b); ++n) { + Real &v0n = v0(b, n, k, j, i); + Real &v1n = v1(b, n, k, j, i); + v0n = g0 * v0n + g1 * v1n; + if constexpr (include_divf) { + v0n += bdt_vol * ((ax1[0] * v0.flux(b, d1, n, k, j, i) - + ax1[1] * v0.flux(b, d1, n, k, j, i + 1)) + + (ax2[0] * v0.flux(b, d2, n, k, j, i) - + ax2[1] * v0.flux(b, d2, n, k, j + multi_d, i)) + + (ax3[0] * v0.flux(b, d3, n, k, j, i) - + ax3[1] * v0.flux(b, d3, n, k + three_d, j, i))); + } + } }); return TaskStatus::complete; } From 365a13db9fd5801e9ac8b2f6e57d385fda19e369 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 17 Dec 2025 14:01:35 -0500 Subject: [PATCH 4/8] Eliminate commented out pure MG solver --- inputs/self-gravity/poly.in | 2 +- src/gravity/gravity.cpp | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/inputs/self-gravity/poly.in b/inputs/self-gravity/poly.in index b24133b8..17292b04 100644 --- a/inputs/self-gravity/poly.in +++ b/inputs/self-gravity/poly.in @@ -87,7 +87,7 @@ gamma = 2.0 precondition = true print_per_step = true max_iterations = 100 -residual_tolerance = 1.0e-20 +residual_tolerance = 1.0e-12 smoother = SRJ2 do_FAS = true units_override = true diff --git a/src/gravity/gravity.cpp b/src/gravity/gravity.cpp index d0e1cb14..670df581 100644 --- a/src/gravity/gravity.cpp +++ b/src/gravity/gravity.cpp @@ -212,19 +212,16 @@ std::shared_ptr Initialize(ParameterInput *pin, {Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::FillGhost}); gravity->AddField(mrhs); - // Multigrid + // Solvers using PoissEq = PoissonEquation; PoissEq eq(pin, "poisson"); params.Add("poisson_equation", eq, parthenon::Params::Mutability::Mutable); - std::shared_ptr psolver; using prolongator_t = parthenon::solvers::ProlongationBlockInteriorZeroDirichlet; using preconditioner_t = parthenon::solvers::MGSolver; psolver = std::make_shared>( "base", "phi", "rhs", pin, block_name, PoissEq(pin, block_name)); - // psolver = std::make_shared>( - // "base", "phi", "rhs", pin, block_name, PoissEq(pin, block_name)); params.Add("solver_pointer", psolver); } From 04ad2b3431702785f4a336e056d9254f9448d94d Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Mon, 5 Jan 2026 15:53:31 -0500 Subject: [PATCH 5/8] Try eliminating zero_phi call --- src/gravity/poisson_driver.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gravity/poisson_driver.cpp b/src/gravity/poisson_driver.cpp index f4a5a480..96521412 100644 --- a/src/gravity/poisson_driver.cpp +++ b/src/gravity/poisson_driver.cpp @@ -63,9 +63,9 @@ void SolvePoisson(TaskCollection &tc, Mesh *pmesh) { md, md_rhs); // Set initial solution guess to zero - auto zero_phi = - tl.AddTask(copy_rhs, TF(solvers::utils::SetToZero), md_phi); - auto setup = psolver->AddSetupTasks(tl, zero_phi, i, pmesh); + // auto zero_phi = + // tl.AddTask(copy_rhs, TF(solvers::utils::SetToZero), md_phi); + auto setup = psolver->AddSetupTasks(tl, copy_rhs, i, pmesh); auto solve = psolver->AddTasks(tl, setup, i, pmesh); // Set BCs after solve From eb5a05761f57a963da8e88c2e5674dfdee489fcf Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 6 Jan 2026 10:03:50 -0500 Subject: [PATCH 6/8] Update polytrope input file --- inputs/self-gravity/poly.in | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/inputs/self-gravity/poly.in b/inputs/self-gravity/poly.in index 17292b04..1df97e65 100644 --- a/inputs/self-gravity/poly.in +++ b/inputs/self-gravity/poly.in @@ -25,17 +25,17 @@ variables = gas.prim.density, & gas.prim.pressure, & grav.phi, & grav.rhs -dt = 1.0e-100 # time increment between outputs +dt = 1.0 # time increment between outputs file_type = hst # HDF5 data dump -dt = 1.0e-100 # time increment between outputs +dt = 1.0 # time increment between outputs -nlim = -1 # cycle limit -tlim = 100.0 # time limit (to be reset by problem generator) -integrator = rk1 # time integration algorithm -ncycle_out = 1 # interval for stdout summary info +nlim = -1 # cycle limit +tlim = 100.0 # time limit (to be reset by problem generator) +integrator = rk2 # time integration algorithm +ncycle_out = 1 # interval for stdout summary info # do_coalesced_comms = true @@ -87,7 +87,7 @@ gamma = 2.0 precondition = true print_per_step = true max_iterations = 100 -residual_tolerance = 1.0e-12 +residual_tolerance = 1.0e-6 smoother = SRJ2 do_FAS = true units_override = true @@ -100,4 +100,4 @@ ix3_bc = zero ox3_bc = zero -iprob = 1 +iprob = 2 From 1ff1d4047f5f4a9e46ede9cd9014ec6766c59933 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 6 Jan 2026 10:47:34 -0500 Subject: [PATCH 7/8] Cleanup commented out code --- src/gravity/poisson_driver.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/gravity/poisson_driver.cpp b/src/gravity/poisson_driver.cpp index 96521412..1beb2121 100644 --- a/src/gravity/poisson_driver.cpp +++ b/src/gravity/poisson_driver.cpp @@ -62,9 +62,7 @@ void SolvePoisson(TaskCollection &tc, Mesh *pmesh) { tl.AddTask(copy_rhs, TF(solvers::utils::CopyData>), md, md_rhs); - // Set initial solution guess to zero - // auto zero_phi = - // tl.AddTask(copy_rhs, TF(solvers::utils::SetToZero), md_phi); + // Solve auto setup = psolver->AddSetupTasks(tl, copy_rhs, i, pmesh); auto solve = psolver->AddTasks(tl, setup, i, pmesh); From b29834c2003ac88e84dffedd53af948bdc0680eb Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Wed, 7 Jan 2026 09:59:55 -0500 Subject: [PATCH 8/8] Bump parth --- external/parthenon | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/parthenon b/external/parthenon index 5a8dea63..01c61ebf 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 5a8dea63e41332500059c5ced7aa318ae1549ed0 +Subproject commit 01c61ebffcb388ca349d23ffe3fa692b93dc7b1e