Add Self-Gravity

inputs/self-gravity/poly.in

@@ -0,0 +1,103 @@
+# ========================================================================================
+#  (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
+#  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>
+problem = polytrope      # name of the pgen
+coordinates = cartesian  # coordinate system
+
+<parthenon/job>
+problem_id = polytrope  # problem ID: basename of output filenames
+
+<parthenon/output1>
+file_type = hdf5  # HDF5 data dump
+variables = gas.prim.density, &
+            gas.prim.velocity, &
+            gas.prim.pressure, &
+            grav.phi, &
+            grav.rhs
+dt        = 1.0   # time increment between outputs
+
+<parthenon/output2>
+file_type = hst   # HDF5 data dump
+dt        = 1.0   # time increment between outputs
+
+<parthenon/time>
+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
+
+<parthenon/mesh>
+# 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
+
+<parthenon/meshblock>
+nx1 = 16
+nx2 = 16
+nx3 = 16
+
+<parthenon/refinement0>
+method = magnitude
+comparator = greater_than
+field = gas.prim.density_0
+refine_tol = 0.1
+derefine_tol = 0.1
+
+<physics>
+gas = true
+self_gravity = true
+
+<gas>
+cfl = 0.9
+reconstruct = plm
+riemann = hllc
+gamma = 2.0
+
+<self_gravity>
+precondition = true
+print_per_step = true
+max_iterations = 100
+residual_tolerance = 1.0e-6
+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
+
+<problem>
+iprob = 2

src/CMakeLists.txt

@@ -40,10 +40,10 @@ set (SRC_LIST
 
   gravity/binary_mass.cpp
   gravity/gravity.cpp
+  gravity/gravity.hpp
   gravity/nbody_gravity.hpp
   gravity/point_mass.cpp
   gravity/uniform.cpp
-  gravity/gravity.hpp
 
   nbody/nbody.cpp
   nbody/nbody.hpp
@@ -64,6 +64,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
@@ -85,6 +86,12 @@ set (SRC_LIST
   rotating_frame/rotating_frame_impl.hpp
   rotating_frame/rotating_frame.hpp
 
+  self_gravity/self_gravity.cpp
+  self_gravity/self_gravity.hpp
+  self_gravity/poisson_driver.cpp
+  self_gravity/poisson_equation.hpp
+  self_gravity/self_gravity.cpp
+
   utils/artemis_utils.cpp
   utils/artemis_utils.hpp
   utils/history.hpp

src/artemis.cpp

@@ -24,6 +24,7 @@
 #include "radiation/moments/moments.hpp"
 #include "radiation/radiation.hpp"
 #include "rotating_frame/rotating_frame.hpp"
+#include "self_gravity/self_gravity.hpp"
 #include "utils/artemis_utils.hpp"
 #include "utils/history.hpp"
 #include "utils/units.hpp"
@@ -80,6 +81,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &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 =
@@ -89,6 +91,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   const bool do_gas = pin->GetOrAddBoolean("physics", "gas", true);
   const bool do_dust = pin->GetOrAddBoolean("physics", "dust", false);
   const bool do_gravity = pin->GetOrAddBoolean("physics", "gravity", false);
+  const bool do_self_gravity = pin->GetOrAddBoolean("physics", "self_gravity", false);
   const bool do_nbody = pin->GetOrAddBoolean("physics", "nbody", false);
   const bool do_rotating_frame = pin->GetOrAddBoolean("physics", "rotating_frame", false);
   const bool do_cooling = pin->GetOrAddBoolean("physics", "cooling", false);
@@ -119,6 +122,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   artemis->AddParam("do_gas", do_gas);
   artemis->AddParam("do_dust", do_dust);
   artemis->AddParam("do_gravity", do_gravity);
+  artemis->AddParam("do_self_gravity", do_self_gravity);
   artemis->AddParam("do_nbody", do_nbody);
   artemis->AddParam("do_rotating_frame", do_rotating_frame);
   artemis->AddParam("do_cooling", do_cooling);
@@ -144,6 +148,8 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   // Call package initializers here
   if (do_nbody) packages.Add(NBody::Initialize(pin.get(), constants));
   if (do_gravity) packages.Add(Gravity::Initialize(pin.get(), constants, packages));
+  if (do_self_gravity)
+    packages.Add(SelfGravity::Initialize(pin.get(), constants, packages));
   if (do_gas) packages.Add(Gas::Initialize(pin.get(), units, constants, packages));
   if (do_dust) packages.Add(Dust::Initialize(pin.get(), units));
   if (do_rotating_frame) packages.Add(RotatingFrame::Initialize(pin.get()));

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)

src/artemis_driver.cpp

@@ -30,6 +30,7 @@
 #include "radiation/imc/imc.hpp"
 #include "radiation/moments/moments.hpp"
 #include "rotating_frame/rotating_frame.hpp"
+#include "self_gravity/self_gravity.hpp"
 #include "utils/integrators/artemis_integrator.hpp"
 
 using namespace parthenon::driver::prelude;
@@ -59,6 +60,7 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in
   do_gas = artemis_pkg->template Param<bool>("do_gas");
   do_dust = artemis_pkg->template Param<bool>("do_dust");
   do_gravity = artemis_pkg->template Param<bool>("do_gravity");
+  do_self_gravity = artemis_pkg->template Param<bool>("do_self_gravity");
   do_rotating_frame = artemis_pkg->template Param<bool>("do_rotating_frame");
   do_shear = artemis_pkg->template Param<bool>("do_shear");
   do_cooling = artemis_pkg->template Param<bool>("do_cooling");
@@ -206,6 +208,9 @@ TaskCollection ArtemisDriver<GEOM>::StepTasks() {
     const Real g1 = integrator->gam1[stage - 1];
     const Real bdt = integrator->beta[stage - 1] * integrator->dt;
 
+    // Compute gravitational potential
+    if (do_self_gravity) SelfGravity::SolvePoisson(tc, pmesh);
+
     TaskRegion &tr = tc.AddRegion(num_partitions);
     for (int i = 0; i < num_partitions; i++) {
       auto &tl = tr[i];
@@ -264,11 +269,19 @@ TaskCollection ArtemisDriver<GEOM>::StepTasks() {
                                  Gravity::ExternalGravity<GEOM>, u0.get(), time, bdt);
       }
 
+      // Apply self-gravity source term
+      TaskID self_gravity_src = gas_coord_src | dust_coord_src | gas_diff_src;
+      if (do_self_gravity) {
+        self_gravity_src =
+            tl.AddTask(gas_coord_src | dust_coord_src | gas_diff_src,
+                       SelfGravity::SelfGravity<GEOM>, u0.get(), time, bdt);
+      }
+
       // Apply rotating frame source term
-      TaskID rframe_src = gravity_src;
+      TaskID rframe_src = gravity_src | self_gravity_src;
       if (do_rotating_frame) {
-        rframe_src = tl.AddTask(gravity_src, RotatingFrame::RotatingFrameForce, u0.get(),
-                                time, bdt);
+        rframe_src = tl.AddTask(gravity_src | self_gravity_src,
+                                RotatingFrame::RotatingFrameForce, u0.get(), time, bdt);
       }
 
       // Apply drag source term

src/artemis_driver.hpp

@@ -53,7 +53,7 @@ class ArtemisDriver : public EvolutionDriver {
   IntegratorPtr_t integrator, nbody_integrator, rad_integrator;
   StateDescriptor *artemis_pkg;
   bool do_gas, do_dust, do_moment, do_imc;
-  bool do_gravity, do_nbody, do_rotating_frame, do_shear;
+  bool do_gravity, do_self_gravity, do_nbody, do_rotating_frame, do_shear;
   bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion;
   const bool is_restart;
 };

src/artemis_params.yaml

@@ -74,6 +74,10 @@ physics:
     _type: "bool"
     _default: "false"
     _description: "Turn on gravity"
+  self_gravity:
+    _type: "bool"
+    _default: "false"
+    _description: "Turn on self-gravity"
   rotating_frame:
     _type: "bool"
     _default: "false"

src/gravity/gravity.cpp

@@ -120,7 +120,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin,
   }
 
   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);

src/gravity/gravity.hpp

@@ -113,8 +113,7 @@ TaskStatus NBodyGravityFixed(MeshData<Real> *md, const Real time, const Real dt)
 template <Coordinates GEOM>
 TaskStatus ExternalGravity(MeshData<Real> *md, const Real time, const Real dt);
 
-KOKKOS_INLINE_FUNCTION
-Real quad_ramp(const Real x) { return SQR(x); }
+KOKKOS_INLINE_FUNCTION Real quad_ramp(const Real x) { return SQR(x); }
 
 } // namespace Gravity
 

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
@@ -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<T>(pmb, pin);
   } else if (name == "lw") {
     lw::ProblemGenerator<T>(pmb, pin);
+  } else if (name == "polytrope") {
+    polytrope::ProblemGenerator<T>(pmb, pin);
   } else if (name == "kh") {
     kh::ProblemGenerator<T>(pmb, pin);
   } else if (name == "rt") {