add custom grav flag to VL integrator and remove print statements

cholla-hydro · evaneschneider · Oct 27, 2023 · Sep 13, 2023 · Sep 13, 2023 · Sep 13, 2023
commit 03c3e3d1065353654de27fbbc0ee5a00195a8800
diff --git a/builds/make.type.static_grav b/builds/make.type.static_grav
@@ -29,4 +29,3 @@ DFLAGS    += -DSTATIC_GRAV
 # Can also add -DSLICES and -DPROJECTIONS
 OUTPUT    ?=  -DOUTPUT -DHDF5
 DFLAGS    += $(OUTPUT)
-DN_OUTPUT_COMPLETE=1
diff --git a/src/global/global.cpp b/src/global/global.cpp
@@ -134,11 +134,6 @@ void parse_params(char *param_file, struct parameters *parms, int argc, char **a
     return;
   }
 
-#ifdef STATIC_GRAV
-  // initialize custom gravity flag to zero
-  parms->custom_grav = 0;
-#endif
-
 #ifdef COSMOLOGY
   // Initialize file name as an empty string
   parms->scale_outputs_file[0] = '\0';

diff --git a/src/global/global.h b/src/global/global.h
@@ -199,7 +199,7 @@ struct parameters {
   int out_float32_GasEnergy = 0;
 #endif
 #ifdef STATIC_GRAV
-  int custom_grav;  // flag to set specific static gravity field
+  int custom_grav = 0;  // flag to set specific static gravity field
 #endif
 #ifdef MHD
   int out_float32_magnetic_x = 0;

diff --git a/src/gravity/static_grav.h b/src/gravity/static_grav.h
@@ -18,8 +18,10 @@ inline __device__ void calc_g_1D(int xid, int x_off, int n_ghost, int custom_gra
 {
   Real x_pos, r_disk, r_halo;
   x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
+  //set gravity field according to parameter file input
   switch (custom_grav) {
-    case 1:
+  case 1:
+    //1D NFW halo & Miyamoto-Nagai disk
       // for disk components, calculate polar r
       // r_disk = 0.220970869121;
       // r_disk = 6.85009694274;
@@ -64,12 +66,12 @@ inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_g
   // positions on the grid
   x_pos = (x_off + xid - n_ghost + 0.5) * dx + xbound;
   y_pos = (y_off + yid - n_ghost + 0.5) * dy + ybound;
-  // for Gresho, also need r & phi
+  // for Gresho and disks, also need r & phi
   r   = sqrt(x_pos * x_pos + y_pos * y_pos);
   phi = atan2(y_pos, x_pos);
   switch (custom_grav) {
     case 1:
-      //      printf("gresho\n");
+      // Gresho vortex
       // set acceleration to balance v_phi in Gresho problem
       if (r < 0.2) {
         *gx = -cos(phi) * 25.0 * r;
@@ -83,19 +85,18 @@ inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_g
       }
       break;
     case 2:
-      // printf("rayleigh talor\n");
+      //Rayleigh-Taylor instability
       *gx = 0;
       *gy = -1;
       break;
     case 3:
-      // printf("keplerian\n");
+      // 2D disk in keplerian rotation
       Real M;
       M   = 1 * MSUN_CGS;
       *gx = -cos(phi) * GN * M / (r * r);
       *gy = -sin(phi) * GN * M / (r * r);
       break;
-    case 4:
-      // printf("disk\n");
+  case 4:
       // set gravitational acceleration for Kuzmin disk + NFW halo
       Real a_d, a_h, a, M_vir, M_d, R_vir, R_d, R_s, M_h, c_vir, x;
       M_vir = 1.0e12;         // viral mass of MW in M_sun
@@ -116,7 +117,6 @@ inline __device__ void calc_g_2D(int xid, int yid, int x_off, int y_off, int n_g
       *gy = -sin(phi) * a;
       break;
     default:
-      // printf("default\n");
       *gx = 0;
       *gy = 0;
   }
@@ -139,31 +139,55 @@ inline __device__ void calc_g_3D(int xid, int yid, int zid, int x_off, int y_off
   r_disk = sqrt(x_pos * x_pos + y_pos * y_pos);
   // for halo, calculate spherical r
   r_halo = sqrt(x_pos * x_pos + y_pos * y_pos + z_pos * z_pos);
+  Real a_disk_r, a_disk_z, a_halo, a_halo_r, a_halo_z;
+  Real M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
   switch (custom_grav) {
     case 1:
+      // Milky way disk model
       // set properties of halo and disk (these must match initial conditions)
-      Real a_disk_r, a_disk_z, a_halo, a_halo_r, a_halo_z;
-      Real M_vir, M_d, R_vir, R_d, z_d, R_h, M_h, c_vir, phi_0_h, x;
-      // MW model
+
       M_vir = 1.0e12;     // viral mass of in M_sun
       M_d   = 6.5e10;     // viral mass of in M_sun
       R_d   = 3.5;        // disk scale length in kpc
       z_d   = 3.5 / 5.0;  // disk scale height in kpc
       R_vir = 261.;       // virial radius in kpc
       c_vir = 20.0;       // halo concentration
-      // M82 model
-      // M_vir = 5.0e10; // viral mass of in M_sun
-      // M_d = 1.0e10; // mass of disk in M_sun
-      // R_d = 0.8; // disk scale length in kpc
-      // z_d = 0.15; // disk scale height in kpc
-      // R_vir = R_d/0.015; // viral radius in kpc
-      // c_vir = 10.0; // halo concentration
-
+
       M_h     = M_vir - M_d;    // halo mass in M_sun
       R_h     = R_vir / c_vir;  // halo scale length in kpc
       phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
       x       = r_halo / R_h;
 
+      // calculate acceleration due to NFW halo & Miyamoto-Nagai disk
+      a_halo   = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
+      a_halo_r = a_halo * (r_disk / r_halo);
+      a_halo_z = a_halo * (z_pos / r_halo);
+      a_disk_r = -GN * M_d * r_disk * pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), -1.5);
+      a_disk_z =
+          -GN * M_d * z_pos * (R_d + sqrt(z_pos * z_pos + z_d * z_d)) /
+          (pow(r_disk * r_disk + pow2(R_d + sqrt(z_pos * z_pos + z_d * z_d)), 1.5) * sqrt(z_pos * z_pos + z_d * z_d));
+
+      // total acceleration is the sum of the halo + disk components
+      *gx = (x_pos / r_disk) * (a_disk_r + a_halo_r);
+      *gy = (y_pos / r_disk) * (a_disk_r + a_halo_r);
+      *gz = a_disk_z + a_halo_z;
+      break;
+  case 2:
+    // M82 model
+    // set properties of halo and disk (these must match initial conditions)
+
+       M_vir = 5.0e10; // viral mass of in M_sun
+       M_d = 1.0e10; // mass of disk in M_sun
+       R_d = 0.8; // disk scale length in kpc
+       z_d = 0.15; // disk scale height in kpc
+       R_vir = R_d/0.015; // viral radius in kpc
+       c_vir = 10.0; // halo concentration
+
+       M_h     = M_vir - M_d;    // halo mass in M_sun
+      R_h     = R_vir / c_vir;  // halo scale length in kpc
+      phi_0_h = GN * M_h / (log(1.0 + c_vir) - c_vir / (1.0 + c_vir));
+      x       = r_halo / R_h;
+
       // calculate acceleration due to NFW halo & Miyamoto-Nagai disk
       a_halo   = -phi_0_h * (log(1 + x) - x / (1 + x)) / (r_halo * r_halo);
       a_halo_r = a_halo * (r_disk / r_halo);

diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp
@@ -149,7 +149,6 @@ void Grid3D::Initialize(struct parameters *P)
 
 #ifdef STATIC_GRAV
   H.custom_grav = P->custom_grav;  // Initialize the custom static gravity flag
-  printf("H.custom_grav is %d\n", H.custom_grav);
   if (H.custom_grav == 0) {
     printf("WARNING: No custom gravity field given. Gravity field will be set to zero.\n");
   }
@@ -439,7 +438,7 @@ Real Grid3D::Update_Grid(void)
   {
 #ifdef CUDA
   #ifdef VL
-    VL_Algorithm_1D_CUDA(C.device, H.nx, x_off, H.n_ghost, H.dx, H.xbound, H.dt, H.n_fields);
+    VL_Algorithm_1D_CUDA(C.device, H.nx, x_off, H.n_ghost, H.dx, H.xbound, H.dt, H.n_fields, H.custom_grav);
   #endif  // VL
   #ifdef SIMPLE
     Simple_Algorithm_1D_CUDA(C.device, H.nx, x_off, H.n_ghost, H.dx, H.xbound, H.dt, H.n_fields, H.custom_grav);
@@ -450,7 +449,7 @@ Real Grid3D::Update_Grid(void)
 #ifdef CUDA
   #ifdef VL
     VL_Algorithm_2D_CUDA(C.device, H.nx, H.ny, x_off, y_off, H.n_ghost, H.dx, H.dy, H.xbound, H.ybound, H.dt,
-                         H.n_fields);
+                         H.n_fields, H.custom_grav);
   #endif  // VL
   #ifdef SIMPLE
     Simple_Algorithm_2D_CUDA(C.device, H.nx, H.ny, x_off, y_off, H.n_ghost, H.dx, H.dy, H.xbound, H.ybound, H.dt,
@@ -462,7 +461,7 @@ Real Grid3D::Update_Grid(void)
 #ifdef CUDA
   #ifdef VL
     VL_Algorithm_3D_CUDA(C.device, C.d_Grav_potential, H.nx, H.ny, H.nz, x_off, y_off, z_off, H.n_ghost, H.dx, H.dy,
-                         H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, density_floor, U_floor,
+                         H.dz, H.xbound, H.ybound, H.zbound, H.dt, H.n_fields, H.custom_grav, density_floor, U_floor,
                          C.Grav_potential);
   #endif  // VL
   #ifdef SIMPLE

diff --git a/src/integrators/VL_1D_cuda.cu b/src/integrators/VL_1D_cuda.cu
@@ -29,7 +29,7 @@ __global__ void Update_Conserved_Variables_1D_half(Real *dev_conserved, Real *de
                                                    int n_fields);
 
 void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Real dx, Real xbound, Real dt,
-                          int n_fields)
+                          int n_fields, int custom_grav)
 {
   // Here, *dev_conserved contains the entire
   // set of conserved variables on the grid
@@ -134,7 +134,7 @@ void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Rea
 
   // Step 6: Update the conserved variable array
   hipLaunchKernelGGL(Update_Conserved_Variables_1D, dimGrid, dimBlock, 0, 0, dev_conserved, F_x, n_cells, x_off,
-                     n_ghost, dx, xbound, dt, gama, n_fields);
+                     n_ghost, dx, xbound, dt, gama, n_fields, custom_grav);
   CudaCheckError();
 
     #ifdef DE

diff --git a/src/integrators/VL_1D_cuda.h b/src/integrators/VL_1D_cuda.h
@@ -9,7 +9,7 @@
     #include "../global/global.h"
 
 void VL_Algorithm_1D_CUDA(Real *d_conserved, int nx, int x_off, int n_ghost, Real dx, Real xbound, Real dt,
-                          int n_fields);
+                          int n_fields, int custom_grav);
 
 void Free_Memory_VL_1D();
 

diff --git a/src/integrators/VL_2D_cuda.cu b/src/integrators/VL_2D_cuda.cu
@@ -26,7 +26,7 @@ __global__ void Update_Conserved_Variables_2D_half(Real *dev_conserved, Real *de
                                                    Real dt, Real gamma, int n_fields);
 
 void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_off, int n_ghost, Real dx, Real dy,
-                          Real xbound, Real ybound, Real dt, int n_fields)
+                          Real xbound, Real ybound, Real dt, int n_fields, int custom_grav)
 {
   // Here, *dev_conserved contains the entire
   // set of conserved variables on the grid
@@ -151,7 +151,7 @@ void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_of
 
   // Step 6: Update the conserved variable array
   hipLaunchKernelGGL(Update_Conserved_Variables_2D, dim2dGrid, dim1dBlock, 0, 0, dev_conserved, F_x, F_y, nx, ny, x_off,
-                     y_off, n_ghost, dx, dy, xbound, ybound, dt, gama, n_fields);
+                     y_off, n_ghost, dx, dy, xbound, ybound, dt, gama, n_fields, custom_grav);
   CudaCheckError();
 
     #ifdef DE

diff --git a/src/integrators/VL_2D_cuda.h b/src/integrators/VL_2D_cuda.h
@@ -9,7 +9,7 @@
     #include "../global/global.h"
 
 void VL_Algorithm_2D_CUDA(Real *d_conserved, int nx, int ny, int x_off, int y_off, int n_ghost, Real dx, Real dy,
-                          Real xbound, Real ybound, Real dt, int n_fields);
+                          Real xbound, Real ybound, Real dt, int n_fields, int custom_grav);
 
 void Free_Memory_VL_2D();
 

diff --git a/src/integrators/VL_3D_cuda.cu b/src/integrators/VL_3D_cuda.cu
@@ -37,7 +37,7 @@ __global__ void Update_Conserved_Variables_3D_half(Real *dev_conserved, Real *de
 
 void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                           int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
-                          Real dt, int n_fields, Real density_floor, Real U_floor, Real *host_grav_potential)
+                          Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor, Real *host_grav_potential)
 {
   // Here, *dev_conserved contains the entire
   // set of conserved variables on the grid
@@ -302,7 +302,7 @@ void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int
   // Step 6: Update the conserved variable array
   hipLaunchKernelGGL(Update_Conserved_Variables_3D, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, Q_Lx, Q_Rx, Q_Ly, Q_Ry,
                      Q_Lz, Q_Rz, F_x, F_y, F_z, nx, ny, nz, x_off, y_off, z_off, n_ghost, dx, dy, dz, xbound, ybound,
-                     zbound, dt, gama, n_fields, density_floor, dev_grav_potential);
+                     zbound, dt, gama, n_fields, custom_grav, density_floor, dev_grav_potential);
   CudaCheckError();
 
   #ifdef MHD

diff --git a/src/integrators/VL_3D_cuda.h b/src/integrators/VL_3D_cuda.h
@@ -10,7 +10,7 @@
 
 void VL_Algorithm_3D_CUDA(Real *d_conserved, Real *d_grav_potential, int nx, int ny, int nz, int x_off, int y_off,
                           int z_off, int n_ghost, Real dx, Real dy, Real dz, Real xbound, Real ybound, Real zbound,
-                          Real dt, int n_fields, Real density_floor, Real U_floor, Real *host_grav_potential);
+                          Real dt, int n_fields, int custom_grav, Real density_floor, Real U_floor, Real *host_grav_potential);
 
 void Free_Memory_VL_3D();