diff --git a/builds/make.type.starblast b/builds/make.type.starblast
index 53df92a9c..582f4cc44 100644
--- a/builds/make.type.starblast
+++ b/builds/make.type.starblast
@@ -11,13 +11,13 @@ DFLAGS += -DANALYSIS
 #DFLAGS += -DPARTICLES_KDK
 
 
-DFLAGS += -DGRAVITY
-DFLAGS += -DGRAVITY_GPU
+#DFLAGS += -DGRAVITY
+#DFLAGS += -DGRAVITY_GPU
 # Use both -DSOR and -DPARIS_GALACTIC to run analytic test and compare solutions
 #DFLAGS += -DSOR
-DFLAGS += -DPARIS_GALACTIC
+#DFLAGS += -DPARIS_GALACTIC
 #DFLAGS += -DGRAVITY_ANALYTIC_COMP
-DFLAGS += -DGRAVITY_5_POINTS_GRADIENT
+#DFLAGS += -DGRAVITY_5_POINTS_GRADIENT
 
 #DFLAGS += -DSTATIC_GRAV
 
diff --git a/src/gravity/grav3D.cpp b/src/gravity/grav3D.cpp
index 866663589..8f29c51a1 100644
--- a/src/gravity/grav3D.cpp
+++ b/src/gravity/grav3D.cpp
@@ -16,9 +16,8 @@
 
 Grav3D::Grav3D(void) {}
 
-void Grav3D::Initialize(Real x_min, Real y_min, Real z_min, Real x_max, Real y_max, Real z_max, Real Lx, Real Ly,
-                        Real Lz, int nx, int ny, int nz, int nx_real, int ny_real, int nz_real, Real dx_real,
-                        Real dy_real, Real dz_real, int n_ghost_pot_offset, struct Parameters *P)
+void Grav3D::Initialize(const SpatialDomainProps &spatial_props, Real Lx, Real Ly, Real Lz, int n_ghost_pot_offset,
+                        Parameters *P)
 {
   // Set Box Size
   Lbox_x = Lx;
@@ -26,29 +25,29 @@ void Grav3D::Initialize(Real x_min, Real y_min, Real z_min, Real x_max, Real y_m
   Lbox_z = Lz;
 
   // Set Box Left boundary positions
-  xMin = x_min;
-  yMin = y_min;
-  zMin = z_min;
+  xMin = spatial_props.xMin;
+  yMin = spatial_props.yMin;
+  zMin = spatial_props.zMin;
 
   // Set Box Right boundary positions
-  xMax = x_max;
-  yMax = y_max;
-  zMax = z_max;
+  xMax = spatial_props.xMax;
+  yMax = spatial_props.yMax;
+  zMax = spatial_props.zMax;
 
   // Set uniform ( dx, dy, dz )
-  dx = dx_real;
-  dy = dy_real;
-  dz = dz_real;
+  dx = spatial_props.dx;
+  dy = spatial_props.dy;
+  dz = spatial_props.dz;
 
   // Set Box Total number of cells
-  nx_total = nx;
-  ny_total = ny;
-  nz_total = nz;
+  nx_total = spatial_props.nx_total;
+  ny_total = spatial_props.ny_total;
+  nz_total = spatial_props.nz_total;
 
   // Set Box local domain number of cells
-  nx_local = nx_real;
-  ny_local = ny_real;
-  nz_local = nz_real;
+  nx_local = spatial_props.nx_local;
+  ny_local = spatial_props.ny_local;
+  nz_local = spatial_props.nz_local;
 
   // Local n_cells without ghost cells
   n_cells = nx_local * ny_local * nz_local;
diff --git a/src/gravity/grav3D.h b/src/gravity/grav3D.h
index 69cf8308a..8e71b6983 100644
--- a/src/gravity/grav3D.h
+++ b/src/gravity/grav3D.h
@@ -4,6 +4,7 @@
 #include <stdio.h>
 
 #include "../global/global.h"
+#include "../grid/spatial_domain_props.h"
 
 #ifdef SOR
   #include "../gravity/potential_SOR_3D.h"
@@ -191,9 +192,8 @@ class Grav3D
 
   /*! \fn void Initialize(int nx_in, int ny_in, int nz_in)
    *  \brief Initialize the grid. */
-  void Initialize(Real x_min, Real y_min, Real z_min, Real x_max, Real y_max, Real z_max, Real Lx, Real Ly, Real Lz,
-                  int nx_total, int ny_total, int nz_total, int nx_real, int ny_real, int nz_real, Real dx_real,
-                  Real dy_real, Real dz_real, int n_ghost_pot_offset, struct Parameters *P);
+  void Initialize(const SpatialDomainProps &spatial_props, Real Lx, Real Ly, Real Lz, int n_ghost_pot_offset,
+                  struct Parameters *P);
 
   void AllocateMemory_CPU(void);
   void Initialize_values_CPU();
diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp
index 7f1f9cc28..f8e3abf92 100644
--- a/src/gravity/gravity_functions.cpp
+++ b/src/gravity/gravity_functions.cpp
@@ -358,9 +358,8 @@ static void printDiff(const Real *p, const Real *q, const int nx, const int ny,
 void Grid3D::Initialize_Gravity(struct Parameters *P)
 {
   chprintf("\nInitializing Gravity... \n");
-  Grav.Initialize(H.xblocal, H.yblocal, H.zblocal, H.xblocal_max, H.yblocal_max, H.zblocal_max, H.xdglobal, H.ydglobal,
-                  H.zdglobal, P->nx, P->ny, P->nz, H.nx_real, H.ny_real, H.nz_real, H.dx, H.dy, H.dz,
-                  H.n_ghost_potential_offset, P);
+  SpatialDomainProps spatial_props = SpatialDomainProps::From_Grid3D(*this, P);
+  Grav.Initialize(spatial_props, H.xdglobal, H.ydglobal, H.zdglobal, H.n_ghost_potential_offset, P);
   chprintf("Gravity Successfully Initialized. \n\n");
 
   if (P->bc_potential_type == 1) {
diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h
index 7b6a7ae59..01635fdf0 100644
--- a/src/grid/grid3D.h
+++ b/src/grid/grid3D.h
@@ -845,13 +845,13 @@ class Grid3D
   void Advance_Particles_KDK_Step2_GPU();
   void Set_Particles_Boundary_GPU(int dir, int side);
   void Set_Particles_Density_Boundaries_Periodic_GPU(int direction, int side);
+  int Load_Particles_Density_Boundary_to_Buffer_GPU(int direction, int side, Real *buffer);
+  void Unload_Particles_Density_Boundary_From_Buffer_GPU(int direction, int side, Real *buffer);
   #endif  // PARTICLES_GPU
   #ifdef GRAVITY_GPU
   void Copy_Potential_From_GPU();
   void Copy_Particles_Density_to_GPU();
   void Copy_Particles_Density_GPU();
-  int Load_Particles_Density_Boundary_to_Buffer_GPU(int direction, int side, Real *buffer);
-  void Unload_Particles_Density_Boundary_From_Buffer_GPU(int direction, int side, Real *buffer);
   #endif  // GRAVITY_GPU
 #endif    // PARTICLES
 
diff --git a/src/grid/spatial_domain_props.cpp b/src/grid/spatial_domain_props.cpp
new file mode 100644
index 000000000..132c897ba
--- /dev/null
+++ b/src/grid/spatial_domain_props.cpp
@@ -0,0 +1,64 @@
+#include "../grid/spatial_domain_props.h"
+
+#include "../gravity/grav3D.h"
+#include "../grid/grid3D.h"
+
+SpatialDomainProps SpatialDomainProps::From_Grav3D(Grav3D& grav)
+{
+  SpatialDomainProps out;
+
+  out.nx_local = grav.nx_local;
+  out.ny_local = grav.ny_local;
+  out.nz_local = grav.nz_local;
+
+  out.nx_total = grav.nx_total;
+  out.ny_total = grav.ny_total;
+  out.nz_total = grav.nz_total;
+
+  out.dx = grav.dx;
+  out.dy = grav.dy;
+  out.dz = grav.dz;
+
+  // Left boundaries of the local domain
+  out.xMin = grav.xMin;
+  out.yMin = grav.yMin;
+  out.zMin = grav.zMin;
+
+  // Right boundaries of the local domain
+  out.xMax = grav.xMax;
+  out.yMax = grav.yMax;
+  out.zMax = grav.zMax;
+
+  return out;
+}
+
+SpatialDomainProps SpatialDomainProps::From_Grid3D(Grid3D& grid, struct Parameters* P)
+{
+  SpatialDomainProps out;
+
+  // Set Box Left boundary positions
+  out.xMin = grid.H.xblocal;  // x_min
+  out.yMin = grid.H.yblocal;  // y_min
+  out.zMin = grid.H.zblocal;  // z_min
+
+  // Set Box Right boundary positions
+  out.xMax = grid.H.xblocal_max;  // x_max;
+  out.yMax = grid.H.yblocal_max;  // y_max;
+  out.zMax = grid.H.zblocal_max;  // z_max;
+
+  // Set uniform ( dx, dy, dz )
+  out.dx = grid.H.dx;  // dx_real;
+  out.dy = grid.H.dy;  // dy_real;
+  out.dz = grid.H.dz;  // dz_real;
+
+  // Set Box Total number of cells
+  out.nx_total = P->nx;  // nx;
+  out.ny_total = P->ny;  // ny;
+  out.nz_total = P->nz;  // nz;
+
+  // Set Box local domain number of cells
+  out.nx_local = grid.H.nx_real;  // nx_real;
+  out.ny_local = grid.H.ny_real;  // ny_real;
+  out.nz_local = grid.H.nz_real;  // nz_real;
+  return out;
+};
\ No newline at end of file
diff --git a/src/grid/spatial_domain_props.h b/src/grid/spatial_domain_props.h
new file mode 100644
index 000000000..6527dd4a2
--- /dev/null
+++ b/src/grid/spatial_domain_props.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include "../global/global.h"
+
+class Grav3D;
+class Grid3D;
+struct Parameters;
+
+/* This is a collection of 15 quantities that appear in 3 other locations throughout the codebase.
+ *
+ * The struct primarily exists to simplify the process of copying these values from one place to
+ * another. (But it may make sense to refactor other parts of the code in terms of this object)
+ */
+struct SpatialDomainProps {
+  // number of cells in the local domain
+  int nx_local, ny_local, nz_local;
+
+  // total number of cells in the entire (global) domain
+  int nx_total, ny_total, nz_total;
+
+  // Left boundaries of the local domain
+  Real xMin, yMin, zMin;
+
+  // Right boundaries of the local domain
+  Real xMax, yMax, zMax;
+
+  // cell widths
+  Real dx, dy, dz;
+
+  static SpatialDomainProps From_Grav3D(Grav3D& grav);
+  static SpatialDomainProps From_Grid3D(Grid3D& grid, Parameters* P);
+};
\ No newline at end of file
diff --git a/src/particles/density_CIC.cpp b/src/particles/density_CIC.cpp
index 428a0e864..f576aafa4 100644
--- a/src/particles/density_CIC.cpp
+++ b/src/particles/density_CIC.cpp
@@ -70,7 +70,7 @@ void Grid3D::Copy_Particles_Density()
   Copy_Particles_Density_to_GPU();
     #endif
   Copy_Particles_Density_GPU();
-  #else
+  #elif defined(GRAVITY)
 
     #ifndef PARALLEL_OMP
   Copy_Particles_Density_function(0, Grav.nz_local);
@@ -90,11 +90,12 @@ void Grid3D::Copy_Particles_Density()
   }
     #endif  // PARALLEL_OMP
 
-  #endif  // GRAVITY_GPU
+  #endif  // GRAVITY_GPU and GRAVITY
 }
 
 void Grid3D::Copy_Particles_Density_function(int g_start, int g_end)
 {
+  #ifdef GRAVITY
   int nx_part, ny_part, nz_part, nGHST;
   nGHST   = Particles.G.n_ghost_particles_grid;
   nx_part = Particles.G.nx_local + 2 * nGHST;
@@ -116,6 +117,7 @@ void Grid3D::Copy_Particles_Density_function(int g_start, int g_end)
       }
     }
   }
+  #endif  // GRAVITY
 }
 
 // Clear the density array: density=0
diff --git a/src/particles/density_boundaries_gpu.cu b/src/particles/density_boundaries_gpu.cu
index fd5c4ddca..756d19b6a 100644
--- a/src/particles/density_boundaries_gpu.cu
+++ b/src/particles/density_boundaries_gpu.cu
@@ -1,9 +1,10 @@
-#if defined(PARTICLES_GPU) && defined(GRAVITY_GPU)
+#if defined(PARTICLES_GPU)
 
   #include <iostream>
 
   #include "../grid/grid3D.h"
   #include "../io/io.h"
+  #include "../utils/error_handling.h"
   #include "particles_3D.h"
 
 __global__ void Set_Particles_Density_Boundaries_Periodic_kernel(int direction, int side, int n_i, int n_j, int nx,
@@ -68,6 +69,9 @@ __global__ void Set_Particles_Density_Boundaries_Periodic_kernel(int direction,
 
 void Grid3D::Set_Particles_Density_Boundaries_Periodic_GPU(int direction, int side)
 {
+  #ifndef GRAVITY_GPU
+  CHOLLA_ERROR("This function should not be invoked when compiled without GPU-Gravity");
+  #endif  // GRAVITY_GPU
   int n_ghost, nx_g, ny_g, nz_g, size, ngrid, n_i, n_j;
   n_ghost = Particles.G.n_ghost_particles_grid;
   nx_g    = Particles.G.nx_local + 2 * n_ghost;
@@ -148,6 +152,9 @@ __global__ void Load_Particles_Density_Boundary_to_Buffer_kernel(int direction,
 
 int Grid3D::Load_Particles_Density_Boundary_to_Buffer_GPU(int direction, int side, Real *buffer)
 {
+    #ifndef GRAVITY_GPU
+  CHOLLA_ERROR("This function should not be invoked when compiled without GPU-Gravity");
+    #endif  // GRAVITY_GPU
   int n_ghost, nx_g, ny_g, nz_g, size_buffer, ngrid, n_i, n_j;
   n_ghost = Particles.G.n_ghost_particles_grid;
   nx_g    = Particles.G.nx_local + 2 * n_ghost;
@@ -236,6 +243,9 @@ __global__ void Unload_Particles_Density_Boundary_to_Buffer_kernel(int direction
 
 void Grid3D::Unload_Particles_Density_Boundary_From_Buffer_GPU(int direction, int side, Real *buffer)
 {
+    #ifndef GRAVITY_GPU
+  CHOLLA_ERROR("This function should not be invoked when compiled without GPU-Gravity");
+    #endif  // GRAVITY_GPU
   int n_ghost, nx_g, ny_g, nz_g, size_buffer, ngrid, n_i, n_j;
   n_ghost = Particles.G.n_ghost_particles_grid;
   nx_g    = Particles.G.nx_local + 2 * n_ghost;
@@ -276,4 +286,4 @@ void Grid3D::Unload_Particles_Density_Boundary_From_Buffer_GPU(int direction, in
 
   #endif  // MPI_CHOLLA
 
-#endif  // PARTICLES_GPU & GRAVITY_GPU
+#endif  // PARTICLES_GPU
diff --git a/src/particles/gravity_CIC.cpp b/src/particles/gravity_CIC.cpp
index 495e7cf33..c6bc5368d 100644
--- a/src/particles/gravity_CIC.cpp
+++ b/src/particles/gravity_CIC.cpp
@@ -46,7 +46,7 @@ void Grid3D::Get_Gravity_Field_Particles()
     #endif  // PARALLEL_OMP
   #endif    // PARTICLES_CPU
 
-  #ifdef PARTICLES_GPU
+  #if defined(PARTICLES_GPU) && defined(GRAVITY_GPU)
   Particles.Get_Gravity_Field_Particles_GPU(Grav.F.potential_h);
   #endif
 }
diff --git a/src/particles/gravity_CIC_gpu.cu b/src/particles/gravity_CIC_gpu.cu
index 4711b1a32..62e30555f 100644
--- a/src/particles/gravity_CIC_gpu.cu
+++ b/src/particles/gravity_CIC_gpu.cu
@@ -108,9 +108,10 @@ void Particles3D::Get_Gravity_Field_Particles_GPU_function(int nx_local, int ny_
                                                            Real *gravity_x_dev, Real *gravity_y_dev,
                                                            Real *gravity_z_dev)
 {
-    #ifndef GRAVITY_GPU
+    #ifdef GRAVITY
+      #ifndef GRAVITY_GPU
   Copy_Potential_To_GPU(potential_host, potential_dev, n_cells_potential);
-    #endif
+      #endif
 
   int nx_g, ny_g, nz_g;
   nx_g = nx_local + 2 * N_GHOST_POTENTIAL;
@@ -133,6 +134,7 @@ void Particles3D::Get_Gravity_Field_Particles_GPU_function(int nx_local, int ny_
                      gravity_y_dev, gravity_z_dev, nx_local, ny_local, nz_local, n_ghost_particles_grid,
                      N_GHOST_POTENTIAL, dx, dy, dz);
   GPU_Error_Check();
+    #endif  // GRAVITY
 }
 
 // Get CIC indexes from the particles positions
diff --git a/src/particles/particles_3D.cpp b/src/particles/particles_3D.cpp
index 9f8d207de..04333fcc9 100644
--- a/src/particles/particles_3D.cpp
+++ b/src/particles/particles_3D.cpp
@@ -29,7 +29,13 @@ void Grid3D::Initialize_Particles(struct Parameters *P)
 {
   chprintf("\nInitializing Particles...\n");
 
-  Particles.Initialize(P, Grav, H.xbound, H.ybound, H.zbound, H.xdglobal, H.ydglobal, H.zdglobal);
+  #ifdef GRAVITY
+  SpatialDomainProps spatial_props = SpatialDomainProps::From_Grav3D(Grav);
+  #else
+  SpatialDomainProps spatial_props = SpatialDomainProps::From_Grid3D(*this, P);
+  #endif
+
+  Particles.Initialize(P, spatial_props, H.xbound, H.ybound, H.zbound, H.xdglobal, H.ydglobal, H.zdglobal);
 
   #if defined(PARTICLES_GPU) && defined(GRAVITY_GPU)
   // Set the GPU array for the particles potential equal to the Gravity GPU
@@ -47,8 +53,8 @@ void Grid3D::Initialize_Particles(struct Parameters *P)
   chprintf("Particles Initialized Successfully. \n\n");
 }
 
-void Particles3D::Initialize(struct Parameters *P, Grav3D &Grav, Real xbound, Real ybound, Real zbound, Real xdglobal,
-                             Real ydglobal, Real zdglobal)
+void Particles3D::Initialize(struct Parameters *P, const SpatialDomainProps &spatial_props, Real xbound, Real ybound,
+                             Real zbound, Real xdglobal, Real ydglobal, Real zdglobal)
 {
   // Initialize local and total number of particles to 0
   n_local         = 0;
@@ -108,27 +114,27 @@ void Particles3D::Initialize(struct Parameters *P, Grav3D &Grav, Real xbound, Re
 
   // Initialize Grid Values
   // Local and total number of cells
-  G.nx_local = Grav.nx_local;
-  G.ny_local = Grav.ny_local;
-  G.nz_local = Grav.nz_local;
-  G.nx_total = Grav.nx_total;
-  G.ny_total = Grav.ny_total;
-  G.nz_total = Grav.nz_total;
+  G.nx_local = spatial_props.nx_local;
+  G.ny_local = spatial_props.ny_local;
+  G.nz_local = spatial_props.nz_local;
+  G.nx_total = spatial_props.nx_total;
+  G.ny_total = spatial_props.ny_total;
+  G.nz_total = spatial_props.nz_total;
 
   // Uniform (dx, dy, dz)
-  G.dx = Grav.dx;
-  G.dy = Grav.dy;
-  G.dz = Grav.dz;
+  G.dx = spatial_props.dx;
+  G.dy = spatial_props.dy;
+  G.dz = spatial_props.dz;
 
   // Left boundaries of the local domain
-  G.xMin = Grav.xMin;
-  G.yMin = Grav.yMin;
-  G.zMin = Grav.zMin;
+  G.xMin = spatial_props.xMin;
+  G.yMin = spatial_props.yMin;
+  G.zMin = spatial_props.zMin;
 
   // Right boundaries of the local domain
-  G.xMax = Grav.xMax;
-  G.yMax = Grav.yMax;
-  G.zMax = Grav.zMax;
+  G.xMax = spatial_props.xMax;
+  G.yMax = spatial_props.yMax;
+  G.zMax = spatial_props.zMax;
 
   // Left boundaries of the global domain
   G.domainMin_x = xbound;
@@ -176,8 +182,16 @@ void Particles3D::Initialize(struct Parameters *P, Grav3D &Grav, Real xbound, Re
     #endif
 
   G.size_blocks_array = 1024 * 128;
-  G.n_cells_potential = (G.nx_local + 2 * N_GHOST_POTENTIAL) * (G.ny_local + 2 * N_GHOST_POTENTIAL) *
-                        (G.nz_local + 2 * N_GHOST_POTENTIAL);
+  {
+    #ifdef N_GHOST_POTENTIAL
+    int n_ghost_pot = N_GHOST_POTENTIAL;
+    #else
+    int n_ghost_pot = 0;
+    #endif
+
+    G.n_cells_potential =
+        (G.nx_local + 2 * n_ghost_pot) * (G.ny_local + 2 * n_ghost_pot) * (G.nz_local + 2 * n_ghost_pot);
+  }
 
     #ifdef SINGLE_PARTICLE_MASS
   mass_dev = NULL;  // This array won't be used
diff --git a/src/particles/particles_3D.h b/src/particles/particles_3D.h
index 2ddd146e7..845cb233f 100644
--- a/src/particles/particles_3D.h
+++ b/src/particles/particles_3D.h
@@ -14,8 +14,11 @@
 
     #include "../global/global.h"
     #include "../gravity/grav3D.h"
+    #include "../grid/spatial_domain_props.h"
 
     #ifdef PARTICLES_GPU
+      #include "../utils/gpu.hpp"  // cudaFree
+
       #define TPB_PARTICLES 1024
       // #define PRINT_GPU_MEMORY
       #define PRINT_MAX_MEMORY_USAGE
@@ -218,8 +221,8 @@ class Particles3D
 
   Particles3D(void);
 
-  void Initialize(struct Parameters *P, Grav3D &Grav, Real xbound, Real ybound, Real zbound, Real xdglobal,
-                  Real ydglobal, Real zdglobal);
+  void Initialize(struct Parameters *P, const SpatialDomainProps &spatial_props, Real xbound, Real ybound, Real zbound,
+                  Real xdglobal, Real ydglobal, Real zdglobal);
 
   void Allocate_Particles_Grid_Field_Real(Real **array_dev, int size);
   void Free_GPU_Array_Real(Real *array);
diff --git a/src/particles/particles_boundaries.cpp b/src/particles/particles_boundaries.cpp
index 96e4f110e..4c1fbedc9 100644
--- a/src/particles/particles_boundaries.cpp
+++ b/src/particles/particles_boundaries.cpp
@@ -96,9 +96,11 @@ void Grid3D::Unload_Particles_From_Buffers_BLOCK(int index, int *flags)
   if (H.TRANSFER_HYDRO_BOUNDARIES) {
     return;
   }
+    #ifdef GRAVITY
   if (Grav.TRANSFER_POTENTIAL_BOUNDARIES) {
     return;
   }
+    #endif
 
   if (index == 0) {
     Unload_Particles_from_Buffer_X0(flags);