Refactor and introduce SphericalOverdensity model

docs/sphinx/Reference/Parameters.md

@@ -43,3 +43,10 @@ These parameters should all be specified in the `[gravity]` parameter table.
 :::{include} param/Gravity.md
 :::
 
+## Spherical Overdensity Model
+
+These parameters should all be specified in the `[model.spherical_overdensity]` parameter table.
+They are required if the {par:param}`init` runtime parameter is set to ``"Spherical_Overpressure_3D"``.
+
+:::{include} param/ModelSphericalOverdensity.md
+:::

docs/sphinx/Reference/param/ModelSphericalOverdensity.md

@@ -0,0 +1,21 @@
+:::{par:parameter} model.spherical_overdensity.overdensity
+
+:Summary: *The desired overdensity*
+:Type: {par:typefmt}`floating-point`
+:Default: *None: must be provided*
+
+Within the spherical overdensity, the density is set to this value.
+
+:::
+
+---
+
+:::{par:parameter} model.spherical_overdensity.bkg_density
+
+:Summary: *The desired background density*
+:Type: {par:typefmt}`floating-point`
+:Default: *None: must be provided*
+
+Outside of the spherical overdensity, the density is set to this value.
+
+:::

docs/sphinx/Reference/param/Required.md

@@ -59,6 +59,47 @@ In 1D and 2D problems, this must be set to 1
 
 ---
 
+
+:::{par:parameter} init
+
+:Summary: Name of initial conditions.
+:Type: {par:typefmt}`string`
+:Default: *None*
+
+The value is case-sensitive.
+
+Current options include:
+- ``"Constant"``
+- ``"Sound_Wave"``
+- ``"Square_Wave"``
+- ``"Riemann"``
+- ``"Shu_Osher"``
+- ``"Blast_1D"``
+- ``"KH"``
+- ``"KH_res_ind"``
+- ``"Rayleigh_Taylor"``
+- ``"Gresho"``
+- ``"Implosion_2D"``
+- ``"Noh_2D"``
+- ``"Noh_3D"``
+- ``"Disk_2D"``
+- ``"Disk_3D"``
+- ``"Disk_3D_particles"``
+- ``"Spherical_Overpressure_3D"``
+- ``"Spherical_Overdensity_3D"``
+- ``"Clouds"``
+- ``"Uniform_Grid"``
+- ``"Zeldovich_Pancake"``
+- ``"Chemistry_Test"``
+- ``"Read_Grid"``
+- ``"Read_Grid_Cat"``
+
+See {repository-file}`src/grid/initial_conditions.cpp` for more information about each option.
+Sample input parameter files for many of these problems can be found in the {repository-dir}`examples` directory.
+:::
+
+---
+
 :::{todo}
 
 Port over the remaining parameters

examples/3D/Spherical_Collapse.txt

@@ -33,3 +33,7 @@ zl_bcnd=1
 zu_bcnd=1
 # path to output directory
 outdir=./
+
+[model.spherical_overdensity]
+bkg_density=0.0005
+overdensity=1.0

examples/scripts/sphere.txt

@@ -32,3 +32,7 @@ zu_bcnd=3
 bc_potential_type=0
 # path to output directory
 outdir=.
+
+[model.spherical_overdensity]
+bkg_density=0.0005
+overdensity=1.0

src/gravity/gravity_boundaries.cu

@@ -213,12 +213,14 @@ void Grid3D::Compute_Potential_Isolated_Boundary(int direction, int side, int bc
   auto [pot_boundary, boundary_buf_props] = tmp;
 
   if (bc_potential_type == 0) {
+    const SphericalOverdensity *overdensity_model = models().try_get<SphericalOverdensity>();
+    CHOLLA_ASSERT(overdensity_model != nullptr, "spherical overdensity model wasn't initialized");
     // Point mass potential GM/r
-    const Real r0       = H.sphere_radius;
-    const Real M        = (H.sphere_density - H.sphere_background_density) * 4.0 * M_PI * r0 * r0 * r0 / 3.0;
-    const Real cm_pos_x = H.sphere_center_x;
-    const Real cm_pos_y = H.sphere_center_y;
-    const Real cm_pos_z = H.sphere_center_z;
+    const Real r0 = overdensity_model->radius;
+    const Real M  = (overdensity_model->overdensity - overdensity_model->bkg_density) * 4.0 * M_PI * r0 * r0 * r0 / 3.0;
+    const Real cm_pos_x = overdensity_model->center_xyz[0];
+    const Real cm_pos_y = overdensity_model->center_xyz[1];
+    const Real cm_pos_z = overdensity_model->center_xyz[2];
     const Real Gconst   = Grav.Gconst;
 
     // define a local function that actually computes the potential

src/gravity/gravity_functions.cpp

@@ -473,19 +473,27 @@ void Grid3D::Compute_Gravitational_Potential(struct Parameters *P)
   Copy_Hydro_Density_to_Gravity();
   #endif
 
+  Real dens_avrg, current_a;
   #ifdef COSMOLOGY
   // If using cosmology, set the gravitational constant to the one in the
   // correct units
   const Real Grav_Constant = Cosmo.cosmo_G;
-  const Real current_a     = Cosmo.current_a;
-  const Real dens_avrg     = Cosmo.rho_0_gas;
+  dens_avrg                = Cosmo.rho_0_gas;
+  current_a                = Cosmo.current_a;
   #else
   const Real Grav_Constant = Grav.Gconst;
   // If slowing the Sphere Collapse problem ( bc_potential_type=0 )
-  const Real dens_avrg = (P->bc_potential_type == 0) ? H.sphere_background_density : 0;
-  const Real r0        = H.sphere_radius;
-  // Re-use current_a as the total mass of the sphere
-  const Real current_a = (H.sphere_density - dens_avrg) * 4.0 * M_PI * r0 * r0 * r0 / 3.0;
+  if (P->bc_potential_type == 0) {
+    const SphericalOverdensity *overdensity_model = models().try_get<SphericalOverdensity>();
+    CHOLLA_ASSERT(overdensity_model != nullptr, "spherical overdensity model wasn't initialized");
+    dens_avrg = overdensity_model->bkg_density;
+    Real r0   = overdensity_model->radius;
+    // Re-use current_a as the total mass of the sphere
+    current_a = (overdensity_model->overdensity - dens_avrg) * 4.0 * M_PI * r0 * r0 * r0 / 3.0;
+  } else {
+    dens_avrg = 0.0;
+    current_a = NAN;  // <- never used
+  }
   #endif
 
   if (!Grav.BC_FLAGS_SET) {

src/grid/grid3D.h

@@ -193,14 +193,6 @@ struct Header {
   // Flag to indicate when to transfer the Conserved boundaries
   bool TRANSFER_HYDRO_BOUNDARIES;
 
-  // Parameters For Spherical Colapse Problem
-  Real sphere_density;
-  Real sphere_radius;
-  Real sphere_background_density;
-  Real sphere_center_x;
-  Real sphere_center_y;
-  Real sphere_center_z;
-
   // only meaningful when GRAVITY and GRAVITY_ANALYTIC_COMP are defined
   bool gas_only_use_static_grav;
 
@@ -606,14 +598,12 @@ class Grid3D
       as per the Noh problem in Liska, 2003, or in Stone, 2008. */
   void Noh_Boundary();
 
-  /*! \fn void Spherical_Overpressure_3D()
-   *  \brief Initialize the grid with a 3D spherical overdensity and
-   * overpressue. */
+  /*! \brief Initialize the grid with a 3D spherical overdensity and overpressue. */
   void Spherical_Overpressure_3D();
 
-  /*! \fn void Spherical_Overpressure_3D()
-   *  \brief Initialize the grid with a 3D spherical overdensity for
-   * gravitational collapse */
+  /*! \brief Initialize the grid with a 3D spherical overdensity for gravitational
+   *  collapse
+   */
   void Spherical_Overdensity_3D();
 
   void Clouds();

src/grid/initial_conditions.cpp

@@ -1209,9 +1209,6 @@ void Grid3D::Disk_2D()
   }
 }
 
-/*! \fn void Spherical_Overpressure_3D()
- *  \brief Spherical overdensity and overpressure causing an spherical explosion
- */
 void Grid3D::Spherical_Overpressure_3D()
 {
   int i, j, k, id;
@@ -1260,32 +1257,29 @@ void Grid3D::Spherical_Overpressure_3D()
   }
 }
 
-/*! \fn void Spherical_Overdensity_3D()
- *  \brief Spherical overdensity for gravitational colapse */
 void Grid3D::Spherical_Overdensity_3D()
 {
+  // in the future, we should either:
+  // - factor out the common logic shared with Grid3D::Spherical_Overpressure_3D, OR
+  // - move this logic to the function where we define the SphericalOverdensity model
   int i, j, k, id;
-  Real x_pos, y_pos, z_pos, r, center_x, center_y, center_z;
-  Real density, pressure, overDensity, overPressure, energy, radius, background_density;
+  Real x_pos, y_pos, z_pos, r;
+  Real density, pressure, overPressure, energy;
   Real vx, vy, vz, v2;
-  center_x = 0.5;
-  center_y = 0.5;
-  center_z = 0.5;
-  // overDensity = 1000 * mu * MP / DENSITY_UNIT; // 100 particles per cm^3
-  overDensity  = 1;
+
+  const SphericalOverdensity *overdensity_model = models().try_get<SphericalOverdensity>();
+  CHOLLA_ASSERT(overdensity_model != nullptr, "spherical overdensity model wasn't initialized");
+  Real center_x           = overdensity_model->center_xyz[0];
+  Real center_y           = overdensity_model->center_xyz[1];
+  Real center_z           = overdensity_model->center_xyz[2];
+  Real overDensity        = overdensity_model->overdensity;
+  Real background_density = overdensity_model->bkg_density;
+  Real radius             = overdensity_model->radius;
+
   overPressure = 0;
   vx           = 0;
   vy           = 0;
   vz           = 0;
-  radius       = 0.2;
-  // background_density = mu * MP / DENSITY_UNIT; // 1 particles per cm^3
-  background_density          = 0.0005;
-  H.sphere_density            = overDensity;
-  H.sphere_radius             = radius;
-  H.sphere_background_density = background_density;
-  H.sphere_center_x           = center_x;
-  H.sphere_center_y           = center_y;
-  H.sphere_center_z           = center_z;
 
   // set the initial values of the conserved variables
   for (k = H.n_ghost; k < H.nz - H.n_ghost; k++) {

src/model/model_collection.cu

@@ -18,6 +18,10 @@ ModelCollection::ModelCollection(ParameterMap& pmap)
   // <model-subtable> will be replaced with the name of the corresponding parameter
   // file subtable
 
+  if (pmap.Contains_Table("model.spherical_overdensity")) {
+    vec_.emplace_back(SphericalOverdensity(pmap));
+  }
+
   // the galaxy_model is still a special case, we will start treating it like a normal
   // case soon
   ClusteredDiskGalaxy galaxy_model = galaxies::make_MW_model();

src/model/model_collection.h

@@ -11,6 +11,7 @@
 #include "../io/ParameterMap.h"
 #include "../utils/error_handling.h"  // always_false
 #include "galaxy/disk_galaxy.h"
+#include "sphere_overdensity/model.h"
 
 /*! \defgroup modelgrp Model Group
  *
@@ -186,15 +187,10 @@
 namespace model_detail
 {
 
-// this is just a placeholder model type until we have 2 or more models
-struct DummyModel {
-  explicit DummyModel(ParameterMap& pmap) {}
-};
-
 // a type-safe union that can represent all model types
 // -> to add a new kind of model, append it to the list of template arguments
 // -> todo: consolidate DiskGalaxy and ClusteredDiskGalaxy into a single class
-using model_variant = std::variant<DummyModel, ClusteredDiskGalaxy>;
+using model_variant = std::variant<ClusteredDiskGalaxy, SphericalOverdensity>;
 
 // define logic to check if a type T is an allowed type of a std::variant
 template <typename T, typename variant>

src/model/sphere_overdensity/model.cpp

@@ -0,0 +1,24 @@
+/*! \file
+ *  \brief Define logic pertaining to the SphericalOverdensity model
+ */
+
+#include "model.h"
+
+#include "../../io/ParameterMap.h"
+
+SphericalOverdensity::SphericalOverdensity(ParameterMap& pmap)
+    // we can consider configuring the following option from pmap at some point in the
+    // the future
+    : radius{0.2}, center_xyz{0.5, 0.5, 0.5}
+{
+  // the following assignments were commented out in the location where we originally
+  // took the initial values from:
+
+  // bkg_density = mu * MP / DENSITY_UNIT; // 1 particles per cm^3)
+  // bkg_density = 0.0005;
+  bkg_density = pmap.value<double>("model.spherical_overdensity.bkg_density");
+
+  // overdensity = 1000 * mu * MP / DENSITY_UNIT; // 100 particles per cm^3
+  // overdensity = 1.0;
+  overdensity = pmap.value<double>("model.spherical_overdensity.overdensity");
+}

src/model/sphere_overdensity/model.h

@@ -0,0 +1,31 @@
+/*! \file
+ *  \brief Define/declare machinery pertaining to spherical overdensity test problem
+ */
+
+#pragma once
+
+#include "../../global/global.h"  // Real
+
+// forward declarations to avoid directly (or indirectly) including the ParameterMap
+// header file
+struct ParameterMap;
+
+/*! \brief A centralized for aggregating properties of the model used with the spherical
+ *  overdensity model.
+ *
+ *  This acts as a model because certain properties need to be known at initialization
+ *  and when updating boundary conditions.
+ *
+ *  \note
+ *  Unless you can get the legacy SOR gravity solver to work properly, it seems highly
+ *  unlikely that this logic will work at all
+ */
+struct SphericalOverdensity {
+  Real bkg_density;
+  Real overdensity;
+  Real radius;
+  Real center_xyz[3];
+
+  /*! \brief primary constructor */
+  explicit SphericalOverdensity(ParameterMap& pmap);
+};

src/system_tests/input_files/tGRAVITYSYSTEMSphericalCollapse_CorrectInputExpectCorrectOutput.txt

@@ -33,3 +33,7 @@ zu_bcnd=1
 # path to output directory
 outdir=./
 bc_potential_type=0
+
+[model.spherical_overdensity]
+bkg_density=0.0005
+overdensity=1.0

src/system_tests/input_files/tPARTICLESSYSTEMSphericalCollapse_CorrectInputExpectCorrectOutput.txt

@@ -35,3 +35,7 @@ zl_bcnd=1
 zu_bcnd=1
 # path to output directory
 outdir=./
+
+[model.spherical_overdensity]
+bkg_density=0.0005
+overdensity=1.0