Multipole boundaries

.gitignore

@@ -11,8 +11,16 @@ parameter_file.txt
 *.h5
 *.bin
 out.*
+*.log
 
 ## dropbox file
 .DS_Store
 
 .remote-sync.json
+
+#ctest files
+*.cmake
+CMakeFiles/
+CMakeCache.txt
+Testing/
+tests/poisson_test/Makefile

Makefile

@@ -13,15 +13,19 @@ OBJS := $(subst .c,.o,$(CFILES)) $(subst .cpp,.o,$(CPPFILES)) $(subst .cu,.o,$(G
 
 #To use GPUs, CUDA must be turned on here
 #Optional error checking can also be enabled
-DFLAGS += -DCUDA #-DCUDA_ERROR_CHECK
+DFLAGS += -DCUDA
+#DFLAGS += -DCUDA_ERROR_CHECK
+
+#Profiling flag to profile only the main loop
+#DFLAGS += -DPROFILING
 
 #To use MPI, DFLAGS must include -DMPI_CHOLLA
 DFLAGS += -DMPI_CHOLLA -DBLOCK
 
 #DFLAGS += -DPRECISION=1
 DFLAGS += -DPRECISION=2
 
-DFLAGS += -DH_CORRECTION
+#DFLAGS += -DH_CORRECTION
 
 # Output
 #DFLAGS += -DBINARY
@@ -43,16 +47,17 @@ DFLAGS += -DPPMP
 DFLAGS += -DHLLC
 
 # Integrator
-#DFLAGS += -DVL
+DFLAGS += -DVL
 #DFLAGS += -DCTU
-DFLAGS += -DSIMPLE
+#DFLAGS += -DSIMPLE
 
 # Dual-Energy Formalism
 #DFLAGS += -DDE
 
 # Apply a minimum value to conserved values
 DFLAGS += -DDENSITY_FLOOR
 DFLAGS += -DTEMPERATURE_FLOOR
+DFLAGS += -DPRESSURE_FLOOR
 
 # Allocate GPU memory only once at the first timestep
 #DFLAGS += -DDYNAMIC_GPU_ALLOC
@@ -65,7 +70,7 @@ DFLAGS += -DTEMPERATURE_FLOOR
 #DFLAGS += -DTILED_INITIAL_CONDITIONS
 
 #Average Slow cell when the cell delta_t is very small
-# DFLAGS += -DAVERAGE_SLOW_CELLS
+#DFLAGS += -DAVERAGE_SLOW_CELLS
 
 #Print Initial Statistics
 DFLAGS += -DPRINT_INITIAL_STATS
@@ -97,10 +102,17 @@ OMP_NUM_THREADS ?= 16
 DFLAGS += -DN_OMP_THREADS=$(OMP_NUM_THREADS)
 #DFLAGS += -DPRINT_OMP_DOMAIN
 
-#Stellar simulation
+# Flags related to the tidal simulation
 DFLAGS += -DTIDES
+# Uses relativistic corrections to the orbit and potential. Otherwise exact Newtonian potential is used
+#DFLAGS += -DTIDES_RELATIVISTIC
+# Outputs the black hole potential, which can be used to compute whether any given fluid cell is bound or unbound
+DFLAGS += -DTIDES_OUTPUT_POTENTIAL_BH
+#Prints the center of mass motion at every step
+#DFLAGS += -DOUTPUT_ALWAYS_COM
+
 
-# Test Poisson solver
+# Test Poisson solver with quasispherical distributions
 #DFLAGS += -DPOISSON_TEST
 
 # Cosmology simulation

make_lux_sor.sh

@@ -1,8 +1,9 @@
 #!/bin/bash
 
+module purge
 module load hdf5/1.10.6
 module load openmpi/4.0.1-cuda
-module load cuda10.1/10.1
+module load cuda10.2/10.2.89
 module load gsl/2.6
 module list
 
@@ -11,4 +12,6 @@ export GRAKLE_HOME='/home/brvillas/code/grackle'
 export POISSON_SOLVER='-DSOR'
 export SUFFIX='.sor'
 make clean
-make -j 16
+make -j 40
+
+source ~/.bashrc

src/VL_3D_cuda.cu

@@ -51,6 +51,7 @@ Real VL_Algorithm_3D_CUDA(Real *host_conserved0, Real *host_conserved1, int nx,
     //printf("Subgrid dimensions set: %d %d %d %d %d %d %d %d %d\n", nx_s, ny_s, nz_s, block1_tot, block2_tot, block3_tot, remainder1, remainder2, remainder3);
     //fflush(stdout);
     block_tot = block1_tot*block2_tot*block3_tot;
+    chprintf("Number of blocks: %i\n", block_tot);
     // number of cells in one subgrid block
     BLOCK_VOL = nx_s*ny_s*nz_s;
     // dimensions for the 1D GPU grid
@@ -115,7 +116,7 @@ Real VL_Algorithm_3D_CUDA(Real *host_conserved0, Real *host_conserved1, int nx,
     CudaSafeCall( cudaMalloc((void**)&dev_dt_array, ngrid*sizeof(Real)) );
     #endif 
 
-    #if defined( GRAVITY ) 
+    #if defined( GRAVITY )
     CudaSafeCall( cudaMalloc((void**)&dev_grav_potential, BLOCK_VOL*sizeof(Real)) );
     #else
     dev_grav_potential = NULL;
@@ -181,6 +182,10 @@ Real VL_Algorithm_3D_CUDA(Real *host_conserved0, Real *host_conserved1, int nx,
     hipLaunchKernelGGL(Update_Conserved_Variables_3D_half, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, dev_conserved_half, F_x, F_y, F_z, nx_s, ny_s, nz_s, n_ghost, dx, dy, dz, 0.5*dt, gama, n_fields, density_floor );
     CudaCheckError();
 
+    #ifdef TEMPERATURE_FLOOR
+    hipLaunchKernelGGL(Apply_Temperature_Floor, dim1dGrid, dim1dBlock, 0, 0, dev_conserved_half, nx_s, ny_s, nz_s, n_ghost, n_fields, U_floor );
+    CudaCheckError();
+    #endif //TEMPERATURE_FLOOR
 
     // Step 4: Construct left and right interface values using updated conserved variables
     #ifdef PCM

src/cuda_mpi_routines.cu

@@ -32,8 +32,8 @@ int initialize_cuda_mpi(int myid, int nprocs)
   //double check
   cudaGetDevice(&i_device);
 
-  // printf("In initialize_cuda_mpi: myid = %d, i_device = %d, n_device = %d\n",myid,i_device,n_device);
-  // fflush(stdout);
+  printf("In initialize_cuda_mpi: myid = %d, i_device = %d, n_device = %d\n",myid,i_device,n_device);
+  fflush(stdout);
 
   return 0;
 

src/global.cpp

@@ -8,8 +8,12 @@
 #include <stdlib.h>
 #include <string.h>
 #include <ctype.h>
-#include"global.h"
+#include "global.h"
+#include "io.h"
 
+#ifdef CUDA
+#include"global_cuda.h"
+#endif
 
 /* Global variables */
 Real gama; // Ratio of specific heats
@@ -270,41 +274,48 @@ parms->scale_outputs_file[0] = '\0';
       parms->Rstar = atof(value);
     else if (strcmp(name, "Mbh")==0)
       parms->Mbh = atof(value);
+    else if (strcmp(name, "relaxRate0")==0)
+      parms->relaxRate0 = atof(value);
+    else if (strcmp(name, "relaxRateBkgnd")==0)
+      parms->relaxRateBkgnd = atof(value);
+    else if (strcmp(name, "rhoAmb")==0)
+      parms->rhoAmb = atof(value);
+    else if (strcmp(name, "pAmb")==0)
+      parms->pAmb = atof(value);
     else if (strcmp(name, "polyN")==0)
       parms->polyN = atof(value);
-		else if (strcmp(name, "rprt")==0)
-			parms->rprt = atof(value);
-		else if (strcmp(name, "tRelaxtDyn")==0)
-			parms->tRelaxtDyn = atof(value);
-		else if (strcmp(name, "r0rt") == 0)
-			parms->r0rt = atof(value);
+    else if (strcmp(name, "rprt")==0)
+      parms->rprt = atof(value);
+    else if (strcmp(name, "tRelaxtDyn")==0)
+      parms->tRelaxtDyn = atof(value);
+    else if (strcmp(name, "r0rt") == 0)
+      parms->r0rt = atof(value);
 #endif
-
 #ifdef POISSON_TEST
-		else if (strcmp(name, "c0")==0)
-			parms->c[0] = atof(value);
-		else if (strcmp(name, "c1")==0)
-			parms->c[1] = atof(value);
-		else if (strcmp(name, "c2")==0)
-			parms->c[2] = atof(value);
-		else if (strcmp(name, "c3")==0)
-			parms->c[3] = atof(value);
-		else if (strcmp(name, "c4")==0)
-			parms->c[4] = atof(value);
-		else if (strcmp(name, "c5")==0)
-			parms->c[5] = atof(value);
-		else if (strcmp(name, "d0")==0)
-			parms->d[0] = atof(value);
-		else if (strcmp(name, "d1")==0)
-			parms->d[1] = atof(value);
-		else if (strcmp(name, "d2")==0)
-			parms->d[2] = atof(value);
-		else if (strcmp(name, "d3")==0)
-			parms->d[3] = atof(value);
-		else if (strcmp(name, "d4")==0)
-			parms->d[4] = atof(value);
-		else if (strcmp(name, "d5")==0)
-			parms->d[5] = atof(value);
+    else if (strcmp(name, "c0")==0)
+      parms->c[0] = atof(value);
+    else if (strcmp(name, "c1")==0)
+      parms->c[1] = atof(value);
+    else if (strcmp(name, "c2")==0)
+      parms->c[2] = atof(value);
+    else if (strcmp(name, "c3")==0)
+      parms->c[3] = atof(value);
+    else if (strcmp(name, "c4")==0)
+      parms->c[4] = atof(value);
+    else if (strcmp(name, "c5")==0)
+      parms->c[5] = atof(value);
+    else if (strcmp(name, "d0")==0)
+      parms->d[0] = atoi(value);
+    else if (strcmp(name, "d1")==0)
+      parms->d[1] = atoi(value);
+    else if (strcmp(name, "d2")==0)
+      parms->d[2] = atoi(value);
+    else if (strcmp(name, "d3")==0)
+      parms->d[3] = atoi(value);
+    else if (strcmp(name, "d4")==0)
+      parms->d[4] = atoi(value);
+    else if (strcmp(name, "d5")==0)
+      parms->d[5] = atoi(value);
 #endif//POISSON_TEST
 
 #ifdef SET_MPI_GRID
@@ -324,3 +335,103 @@ parms->scale_outputs_file[0] = '\0';
   /* Close file */
   fclose (fp);
 }
+
+void printHydroParams(){
+
+  chprintf("\nHydro solver parameters:\n");
+
+//Time integrator
+  chprintf(" Integrator: ");
+  #ifdef CTU
+  chprintf("CTU");
+  #elif defined VL
+  chprintf("VL");
+  #elif defined SIMPLE
+  chprintf("Simple");
+  #else
+  chprintf("not recognized");
+  #endif
+
+  chprintf("\n");
+
+//Reconstruction
+  chprintf(" Reconstruction: ");
+  #ifdef PCM
+  chprintf("PCM");
+  #elif defined PLMP
+  chprintf("PLMP");
+  #elif defined PPLMC
+  chprintf("PPLMC");
+  #elif defined PPMP
+  chprintf("PPMP");
+  #elif defined PPMC
+  chprintf("PPMC");
+  #else
+  chprintf("not recognized");
+  #endif
+
+  chprintf("\n");
+
+//Riemann solver
+  chprintf(" Riemann solver: ");
+  #ifdef EXACT
+  chprintf("exact");
+  #elif defined ROE
+  chprintf("Roe");
+  #elif defined HLLC
+  chprintf("HLLC");
+  #else
+  chprintf("not recognized");
+  #endif
+
+  chprintf("\n");
+
+//H correction
+  chprintf(" H correction: ");
+  #ifdef H_CORRECTION
+  chprintf("enabled");
+  #else
+  chprintf("disabled");
+  #endif
+
+  chprintf("\n");
+
+  chprintf(" CFL: %f\n", C_cfl);
+  chprintf(" Floors:\n");
+  chprintf("  T  : %.10e\n", TEMP_FLOOR);
+  chprintf("  rho: %.10e\n", DENS_FLOOR);
+  chprintf("  P  : %.10e\n", PRES_FLOOR);
+
+}
+
+void printMemoryUsageGPU(){
+
+  size_t free_bytes, total_bytes;
+  cudaError_t cuda_status;
+
+  cuda_status = cudaMemGetInfo( &free_bytes, &total_bytes );
+  if ( cudaSuccess != cuda_status ) printf("Error: cudaMemGetInfo failed, %s \n", cudaGetErrorString(cuda_status) );
+
+  double free_db = (double)free_bytes ;
+  double total_db = (double)total_bytes ;
+  double used_db = total_db - free_db ;
+/*
+  #ifdef MPI_CHOLLA
+  MPI_Allreduce(MPI_IN_PLACE, &used_db, 1, MPI_CHREAL, MPI_MAX, world);
+  #endif
+*/
+
+  #ifdef MPI_CHOLLA
+  char name[MPI_MAX_PROCESSOR_NAME];
+  int len, i_device, n_device;
+  MPI_Get_processor_name( name, &len );
+  cudaGetDeviceCount(&n_device);
+  cudaGetDevice(&i_device);
+
+  printf("Node %s, GPU %i/%d memory usage: %f/%f MB\n", name, i_device, n_device, used_db/1024.0/1024.0, total_db/1024.0/1024.0);
+  MPI_Barrier(MPI_COMM_WORLD);
+  #else
+  printf("Memory usage: %f/%f MB\n", used_db/1024./1024., total_db/1024./1024.);
+  #endif
+
+}

src/global.h

@@ -50,8 +50,9 @@ typedef double Real;
 #define LOG_FILE_NAME "run_output.log"
 
 //Conserved Floor Values
-#define TEMP_FLOOR 1e0
-#define DENS_FLOOR 1e-18
+#define TEMP_FLOOR 0.
+#define DENS_FLOOR (1.e-15)
+#define PRES_FLOOR (1.e-3)
 
 //Parameter for Enzo dual Energy Condition
 #define DE_ETA_1 0.001 //Ratio of U to E for wich  Inetrnal Energy is used to compute the Pressure
@@ -236,15 +237,19 @@ struct parameters
   Real Mstar;
   Real Rstar;
   Real Mbh;
-	Real tRelaxtDyn;
+  Real pAmb;
+  Real rhoAmb;
+  Real tRelaxtDyn;
   Real polyN;
-	Real rprt;
-	Real r0rt;
-#endif//TIDES
+  Real rprt;
+  Real r0rt;
+  Real relaxRate0;
+  Real relaxRateBkgnd;
+#endif
 
 #ifdef POISSON_TEST
-	Real c[6];
-	int d[6];
+  Real c[6];
+  int d[6];
 #endif
 
 #ifdef COSMOLOGY
@@ -272,5 +277,8 @@ struct parameters
  *  \brief Reads the parameters in the given file into a structure. */
 extern void parse_params (char *param_file, struct parameters * parms);
 
+extern void printHydroParams();
+extern void printMemoryUsageGPU();
+
 
 #endif //GLOBAL_H