diff --git a/.gitignore b/.gitignore index 0325338be..ed7e80934 100644 --- a/.gitignore +++ b/.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 diff --git a/Makefile b/Makefile index c489e4652..d76eed754 100644 --- a/Makefile +++ b/Makefile @@ -13,7 +13,11 @@ 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 @@ -21,7 +25,7 @@ DFLAGS += -DMPI_CHOLLA -DBLOCK #DFLAGS += -DPRECISION=1 DFLAGS += -DPRECISION=2 -DFLAGS += -DH_CORRECTION +#DFLAGS += -DH_CORRECTION # Output #DFLAGS += -DBINARY @@ -43,9 +47,9 @@ DFLAGS += -DPPMP DFLAGS += -DHLLC # Integrator -#DFLAGS += -DVL +DFLAGS += -DVL #DFLAGS += -DCTU -DFLAGS += -DSIMPLE +#DFLAGS += -DSIMPLE # Dual-Energy Formalism #DFLAGS += -DDE @@ -53,6 +57,7 @@ DFLAGS += -DSIMPLE # 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 diff --git a/make_lux_sor.sh b/make_lux_sor.sh index 1be033323..78335bd1f 100644 --- a/make_lux_sor.sh +++ b/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 diff --git a/src/VL_3D_cuda.cu b/src/VL_3D_cuda.cu index 3d5ef41e0..1fbb47973 100644 --- a/src/VL_3D_cuda.cu +++ b/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 diff --git a/src/cuda_mpi_routines.cu b/src/cuda_mpi_routines.cu index fd00dbb5d..4dfcc3642 100644 --- a/src/cuda_mpi_routines.cu +++ b/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; diff --git a/src/global.cpp b/src/global.cpp index 85c32632e..9b386618c 100644 --- a/src/global.cpp +++ b/src/global.cpp @@ -8,8 +8,12 @@ #include #include #include -#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 + +} diff --git a/src/global.h b/src/global.h index 20e130909..cdb452383 100644 --- a/src/global.h +++ b/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 diff --git a/src/gravity/grav3D.cpp b/src/gravity/grav3D.cpp index fca6e1265..d360813c0 100644 --- a/src/gravity/grav3D.cpp +++ b/src/gravity/grav3D.cpp @@ -13,11 +13,9 @@ #include "../parallel_omp.h" #endif - - Grav3D::Grav3D( void ){} -void Grav3D::Initialize( Real x_min, Real y_min, Real z_min, 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( Real x_min, Real y_min, Real z_min, Real Lx, Real Ly, Real Lz, int nx, int ny, int nz, int nx_real, int ny_real, int nz_real, int n_ghost, Real dx_real, Real dy_real, Real dz_real, int n_ghost_pot_offset, struct Header H, struct parameters *P ) { //Set Box Size @@ -51,7 +49,7 @@ void Grav3D::Initialize( Real x_min, Real y_min, Real z_min, Real Lx, Real Ly, R n_cells_potential = ( nx_local + 2*N_GHOST_POTENTIAL ) * ( ny_local + 2*N_GHOST_POTENTIAL ) * ( nz_local + 2*N_GHOST_POTENTIAL ); //Set Initial and dt used for the extrapolation of the potential; - //The first timestep the potetential in not extrapolated ( INITIAL = TRUE ) + //The first timestep the potential in not extrapolated ( INITIAL = TRUE ) INITIAL = true; dt_prev = 0; dt_now = 0; @@ -67,16 +65,21 @@ void Grav3D::Initialize( Real x_min, Real y_min, Real z_min, Real Lx, Real Ly, R //Set the Gravitational Constant ( units must be consistent ) Gconst = GN; - if (strcmp(P->init, "Spherical_Overdensity_3D")==0 || strcmp(P->init, "poissonTest") == 0){ + if (strcmp(P->init, "Spherical_Overdensity_3D")==0){ Gconst = 1; - chprintf("WARNING: Using Gravitational Constant G=1.\n"); +// chprintf("WARNING: Using Gravitational Constant G=1.\n"); } + + #ifdef POISSON_TEST + Gconst = 1; + #endif #ifdef TIDES Gconst = G_CGS; - chprintf("WARNING: Using Gravitational Constant in cgs units.\n"); - #endif//TIDES - + #endif + + chprintf(" Using G = %.10e\n", Gconst); + //Flag to transfer the Potential boundaries TRANSFER_POTENTIAL_BOUNDARIES = false; @@ -87,10 +90,11 @@ void Grav3D::Initialize( Real x_min, Real y_min, Real z_min, Real Lx, Real Ly, R Initialize_values_CPU(); - chprintf( "Gravity Initialized: \n Lbox: %0.2f %0.2f %0.2f \n Local: %d %d %d \n Global: %d %d %d \n", - Lbox_x, Lbox_y, Lbox_z, nx_local, ny_local, nz_local, nx_total, ny_total, nz_total ); +// chprintf( "Gravity Initialized: \n Lbox: %0.2f %0.2f %0.2f \n Local: %d %d %d \n Global: %d %d %d \n", +// Lbox_x, Lbox_y, Lbox_z, nx_local, ny_local, nz_local, nx_total, ny_total, nz_total ); +// chprintf("Gravity initialized.\n"); - chprintf( " dx:%f dy:%f dz:%f\n", dx, dy, dz ); +// chprintf( " dx:%f dy:%f dz:%f\n", dx, dy, dz ); chprintf( " N ghost potential: %d\n", N_GHOST_POTENTIAL); chprintf( " N ghost offset: %d\n", n_ghost_pot_offset); @@ -127,6 +131,19 @@ void Grav3D::AllocateMemory_CPU(void) F.pot_boundary_z0 = (Real *) malloc(N_GHOST_POTENTIAL*nx_local*ny_local*sizeof(Real)); //array for the potential isolated boundary F.pot_boundary_z1 = (Real *) malloc(N_GHOST_POTENTIAL*nx_local*ny_local*sizeof(Real)); #endif + + #if defined TIDES || defined POISSON_TEST +//Real and imaginary parts of the multipole moments of the density distribution + ReQ = (Real *) malloc( sizeof(Real) * (LMAX + 1) * ( LMAX + 2 ) / 2); + ImQ = (Real *) malloc( sizeof(Real) * (LMAX + 1) * ( LMAX + 2 ) / 2); + Qblocks = ceil( n_cells / QTPB ); + centerBlocks = ceil( ( nx_local * ny_local * nz_local ) / CENTERTPB ); + bufferReQ = (Real *) malloc( sizeof(Real) * Qblocks * (LMAX + 1) * ( LMAX + 2 ) / 2 ); + bufferImQ = (Real *) malloc( sizeof(Real) * Qblocks * (LMAX + 1) * ( LMAX + 2 ) / 2 ); + bufferCenter = (Real *) malloc( sizeof(Real) * centerBlocks * 3 ); + bufferTotrhosq = (Real *) malloc( sizeof(Real) * centerBlocks ); + #endif + } void Grav3D::Set_Boundary_Flags( int *flags ){ @@ -144,6 +161,12 @@ void Grav3D::Initialize_values_CPU(void){ F.potential_h[id_pot] = 0; F.potential_1_h[id_pot] = 0; } + + for ( int i = 0; i < ( 1 + LMAX ) * ( 2 + LMAX ) / 2; i++ ){ + ReQ[i] = 0.; + ImQ[i] = 0.; + } + } void Grav3D::FreeMemory_CPU(void) diff --git a/src/gravity/grav3D.h b/src/gravity/grav3D.h index cabecaacf..a86a36855 100644 --- a/src/gravity/grav3D.h +++ b/src/gravity/grav3D.h @@ -4,6 +4,12 @@ #include #include"../global.h" +#if defined TIDES || defined POISSON_TEST +#define LMAX (5) +#define QTPB (128) +#define CENTERTPB (1024) +#endif + #ifdef PFFT #include"potential_PFFT_3D.h" #endif @@ -14,6 +20,7 @@ #ifdef SOR #include"potential_SOR_3D.h" +#define SOREPSILON (1.e-8) #endif #ifdef PARIS @@ -107,6 +114,24 @@ class Grav3D #ifdef SOR Potential_SOR_3D Poisson_solver; + + #if defined TIDES || defined POISSON_TEST + Real *ReQ; + Real *ImQ; + Real *bufferReQ; + Real *bufferImQ; + Real *bufferCenter; + Real *bufferTotrhosq; + int Qblocks, centerBlocks; + Real center[3]; + int Qidx(int cidx, int l, int m); + void fillLegP(Real* legP, Real x); + +//TODO: Figure out a better way to pass the density so that we don't have to copy it both for gravity and for the hydro. Right now we're wasting space but it's fine. + Real *dev_rho, *dev_center, *dev_bounds, *dev_dx, *dev_partialReQ, *dev_partialImQ, *dev_partialCenter, *dev_partialTotrhosq; + int *dev_n; + #endif + #endif #ifdef PARIS @@ -154,7 +179,7 @@ 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 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( Real x_min, Real y_min, Real z_min, Real Lx, Real Ly, Real Lz, int nx_total, int ny_total, int nz_total, int nx_real, int ny_real, int nz_real, int n_ghost, Real dx_real, Real dy_real, Real dz_real, int n_ghost_pot_offset, struct Header H, struct parameters *P); void AllocateMemory_CPU(void); void Initialize_values_CPU(); @@ -170,6 +195,12 @@ class Grav3D void Copy_Isolated_Boundaries_To_GPU( struct parameters *P ); #endif + #if defined POISSON_TEST || defined TIDES + void AllocateMemoryBoundaries_GPU(); + void CopyDomainPropertiesToGPU(Real *bounds_local, int *n_local_real, Real *dxi); + void FreeMemoryBoundaries_GPU(); + #endif + }; diff --git a/src/gravity/gravity_boundaries.cpp b/src/gravity/gravity_boundaries.cpp index 322608aed..3b949bb41 100644 --- a/src/gravity/gravity_boundaries.cpp +++ b/src/gravity/gravity_boundaries.cpp @@ -76,7 +76,7 @@ void Grid3D::Set_Potential_Boundaries_Isolated( int direction, int side, int *fl if ( side == 0 ) id_grid = (i+nGHST) + (k)*nx_g + (j+nGHST)*nx_g*ny_g; if ( side == 1 ) id_grid = (i+nGHST) + (k+ny_local+nGHST)*nx_g + (j+nGHST)*nx_g*ny_g; } - if ( direction == 1 ){ + if ( direction == 2 ){ if ( side == 0 ) id_grid = (i+nGHST) + (j+nGHST)*nx_g + (k)*nx_g*ny_g; if ( side == 1 ) id_grid = (i+nGHST) + (j+nGHST)*nx_g + (k+nz_local+nGHST)*nx_g*ny_g; } @@ -100,8 +100,6 @@ void Grid3D::Compute_Potential_Isolated_Boundary( int direction, int side, int Ly_local = Grav.ny_local * Grav.dy; Lz_local = Grav.nz_local * Grav.dz; - - #ifdef GRAV_ISOLATED_BOUNDARY_X if ( direction == 0 ){ domain_l = Grav.xMin; @@ -129,22 +127,12 @@ void Grid3D::Compute_Potential_Isolated_Boundary( int direction, int side, int if ( side == 1 ) pot_boundary = Grav.F.pot_boundary_z1; } #endif - - Real M, cm_pos_x, cm_pos_y, cm_pos_z, pos_x, pos_y, pos_z, r, delta_x, delta_y, delta_z; - M = 0.1005; - cm_pos_x = 0.; - cm_pos_y = 0.; - cm_pos_z = 0.; - - #ifdef TIDES - M = S.Mstar; - cm_pos_x = 0.0; - cm_pos_y = 0.0; - cm_pos_z = 0.0; - #endif//TIDES - + int i, j, k, id; - Real pot_val; + Real pos[3], r, pot_val; + #if defined POISSON_TEST || defined TIDES + Real phi, theta, Ylmfac, lfac; + #endif for ( k=0; knx, 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 ); + Grav.Initialize( H.xblocal, H.yblocal, H.zblocal, H.xdglobal, H.ydglobal, H.zdglobal, P->nx, P->ny, P->nz, H.nx_real, H.ny_real, H.nz_real, H.n_ghost, H.dx, H.dy, H.dz, H.n_ghost_potential_offset, H, P ); chprintf( "Gravity Successfully Initialized. \n\n"); #ifdef PARIS_TEST @@ -353,7 +354,22 @@ void Grid3D::Compute_Gravitational_Potential( struct parameters *P ){ // #endif Grav.BC_FLAGS_SET = true; } - + + +//If doing SOR, copy the density early. The reason is that if doing isolated boundaries with the multipole expansion, we need to know the entire density field to compute the boundaries. Previously the boundaries were computed independently of the rest of the solution. The copy of the density field occurred inside Get_Potential_SOR + #ifdef SOR + Grav.Poisson_solver.Copy_Density_To_GPU(Grav.Poisson_solver.n_cells_local, Grav.F.density_h); + setMoments(); + Grav.Poisson_solver.Convert_Density_To_RHS(Grav_Constant); + if ( !Grav.Poisson_solver.potential_initialized ){ + chprintf( "SOR: Initializing Potential \n"); + Grav.Poisson_solver.Initialize_Potential( Grav.Poisson_solver.nx_local, Grav.Poisson_solver.ny_local, Grav.Poisson_solver.nz_local, Grav.Poisson_solver.n_ghost, Grav.Poisson_solver.F.potential_d, Grav.Poisson_solver.F.density_d ); + Grav.Poisson_solver.potential_initialized = true; + } + +// Grav.Poisson_solver.Copy_Input_And_Initialize( Grav.F.density_h, Grav_Constant, dens_avrg, current_a ); + #endif + #ifdef GRAV_ISOLATED_BOUNDARY_X if ( Grav.boundary_flags[0] == 3 ) Compute_Potential_Boundaries_Isolated(0); if ( Grav.boundary_flags[1] == 3 ) Compute_Potential_Boundaries_Isolated(1); @@ -370,7 +386,6 @@ void Grid3D::Compute_Gravitational_Potential( struct parameters *P ){ // chprintf("Isolated Z\n"); #endif - //Solve Poisson Equation to compute the potential //Poisson Equation: laplacian( phi ) = 4 * pi * G / scale_factor * ( dens - dens_average ) #ifdef SOR @@ -481,7 +496,11 @@ void Grid3D::Extrapolate_Grav_Potential_Function( int g_start, int g_end ){ int nGHST = n_ghost_grid - N_GHOST_POTENTIAL; Real pot_now, pot_prev, pot_extrp; int k, j, i, id_pot, id_grid; - Real posx, posy, posz; + + #ifdef TIDES + Real x[3]; + #endif + for ( k=g_start; k l || l > LMAX || m < 0 ){ + printf("Wrong parameters!\n"); + return -1; + } + + int stride = ( 1 + LMAX ) * ( 2 + LMAX ) / 2; + int substride = l * ( l + 1 ) / 2; + return cidx * stride + substride + m; + +} + +//Same but for the device +__device__ int dQidx(int cidx, int l, int m){ + + if ( m > l || l > LMAX || m < 0 ){ + printf("Wrong parameters!\n"); + return -1; + } + + int stride = ( 1 + LMAX ) * ( 2 + LMAX ) / 2; + int substride = l * ( l + 1 ) / 2; + return cidx * stride + substride + m; + +} + +//Recursively computes Legendre polynomials up to order LMAX. Returns 1D array +__device__ void fillLegP(Real* legP, Real x) +{ + + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + legP[dQidx(0,l,m)] = 0.; + } + } + +//Initial polynomial for recursion relations + legP[dQidx(0,0,0)] = 1./sqrt(4.*M_PI); + +//Diagonal + for( int m = 1; m <= LMAX; m++) + { + legP[dQidx(0,m,m)] = - sqrt( 1. + 1. / 2. / m ) * sqrt( 1. - x * x ) * legP[dQidx(0,m-1,m-1)]; + } + + for( int m = 0; m < LMAX; m++) + { + legP[dQidx(0,m+1,m)] = sqrt( 2. * m + 3. ) * x * legP[dQidx(0,m,m)]; + } + + for( int m = 0; m <= LMAX; m++){ + for( int l = m + 2; l <= LMAX; l++){ + Real c1 = sqrt( ((2.0*l+1)*(2.0*l-1)) / ((l+m)*(l-m))); + Real c2 = sqrt( (2.0*l+1)*(l-m-1.0)*(l+m-1.0) / ((2.0*l-3)*(l-m)*(l+m))); + legP[dQidx(0,l,m)] = c1 * x * legP[dQidx(0,l-1,m)] - c2 * legP[dQidx(0,l-2,m)]; + } + } + +} + +void Grav3D::fillLegP(Real* legP, Real x){ + + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + legP[Qidx(0,l,m)] = 0.; + } + } + +//Initial polynomial for recursion relations + legP[Qidx(0,0,0)] = 1./sqrt(4.*M_PI); + +//Diagonal + for( int m = 1; m <= LMAX; m++) + { + legP[Qidx(0,m,m)] = - sqrt( 1. + 1. / 2. / m ) * sqrt( 1. - x * x ) * legP[Qidx(0,m-1,m-1)]; + } + + for( int m = 0; m < LMAX; m++) + { + legP[Qidx(0,m+1,m)] = sqrt( 2. * m + 3. ) * x * legP[Qidx(0,m,m)]; + } + + for( int m = 0; m <= LMAX; m++){ + for( int l = m + 2; l <= LMAX; l++){ + Real c1 = sqrt( ((2.0*l+1)*(2.0*l-1)) / ((l+m)*(l-m))); + Real c2 = sqrt( (2.0*l+1)*(l-m-1.0)*(l+m-1.0) / ((2.0*l-3)*(l-m)*(l+m))); + legP[Qidx(0,l,m)] = c1 * x * legP[Qidx(0,l-1,m)] - c2 * legP[Qidx(0,l-2,m)]; + } + } + +} + +__device__ int tidFake(int tid_x, int tid_y, int tid_z, int n_ghost, int *n){ + + int tid = ( tid_z + n_ghost ) * n[0] * n[1] + ( tid_y + n_ghost ) * n[0] + ( tid_x + n_ghost ); + #ifdef POISSON_TEST + int tid_z_fake = tid / ( n[0] * n[1]); + int tid_y_fake = ( tid - tid_z_fake * n[0] * n[1] ) / n[0]; + int tid_x_fake = tid - tid_z_fake * n[0] * n[1] - tid_y_fake * n[0]; + if ( tid_x_fake < n_ghost || tid_y_fake < n_ghost || tid_z_fake < n_ghost || tid_z_fake > n[2] - n_ghost || tid_y_fake > n[1] - n_ghost || tid_x_fake > n[0] - n_ghost || tid >= n[0] * n[1] * n[2] || tid_z_fake != tid_z + n_ghost || tid_y_fake != tid_y + n_ghost || tid_x_fake != tid_x + n_ghost){ + printf("Something wrong in cell mapping.\n"); + } + #endif + return tid; +} + +__global__ void QlmKernel(Real *rho, Real *center, Real *bounds, Real *dx, Real rmpole, int *n, int n_ghost, Real *partialReQ, Real *partialImQ){ + + __shared__ Real ReQ[QTPB * (1 + LMAX ) * (2 + LMAX ) / 2]; + __shared__ Real ImQ[QTPB * (1 + LMAX ) * (2 + LMAX ) / 2]; + + for ( int i = threadIdx.x * ( 1 + LMAX ) * ( 2 + LMAX ) / 2; i < ( threadIdx.x + 1 ) * ( 1 + LMAX ) * ( 2 + LMAX ) / 2; i++ ){ + ReQ[i] = 0.; + ImQ[i] = 0.; + } + + int nreal[3]; + for ( int i = 0; i < 3; i++ ) nreal[i] = n[i] - 2 * n_ghost; + int nrealcells = nreal[0] * nreal[1] * nreal[2]; + + int tid = threadIdx.x + blockIdx.x * blockDim.x; + int tid_z = tid / ( nreal[0] * nreal[1] ); + int tid_y = ( tid - tid_z * nreal[0] * nreal[1] ) / nreal[0]; + int tid_x = tid - tid_z * nreal[0] * nreal[1] - tid_y * nreal[0]; + int cidx = threadIdx.x; + int stride = blockDim.x * gridDim.x; + + Real r, phi, fac, pos[3], dev_legP[(1+LMAX)*(2+LMAX)/2]; + + while ( tid < nrealcells ){ + + tid_z = tid / ( nreal[0] * nreal[1] ); + tid_y = ( tid - tid_z * nreal[0] * nreal[1] ) / nreal[0]; + tid_x = tid - tid_z * nreal[0] * nreal[1] - tid_y * nreal[0]; + + pos[0] = bounds[0] + dx[0] * ( tid_x + 0.5) - center[0]; + pos[1] = bounds[1] + dx[1] * ( tid_y + 0.5) - center[1]; + pos[2] = bounds[2] + dx[2] * ( tid_z + 0.5) - center[2]; + r = sqrt( pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2] ); + + if ( r < rmpole ){ + phi = atan2(pos[1], pos[0]); + + fillLegP(dev_legP, pos[2] / r); + + for ( int l = 0; l <= LMAX; l++ ){ + fac = pow(r, l) * rho[tidFake(tid_x, tid_y, tid_z, n_ghost, n)]; + + for ( int m = 0; m <= l; m++ ){ + ReQ[dQidx(cidx, l, m)] += dev_legP[dQidx(0,l,m)] * fac * cos(m * phi); + ImQ[dQidx(cidx, l, m)] += dev_legP[dQidx(0,l,m)] * fac * sin(m * phi); + } + } + } + tid += stride; + } + + __syncthreads(); + + int i = blockDim.x / 2; + while ( i > 0 ){ + if ( cidx < i){ + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + ReQ[dQidx(cidx, l, m)] += ReQ[dQidx(cidx + i, l, m)]; + ImQ[dQidx(cidx, l, m)] += ImQ[dQidx(cidx + i, l, m)]; + } + } + } + __syncthreads(); + i /= 2; + } + + if ( cidx == 0 ){ + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + partialReQ[dQidx(blockIdx.x, l, m)] = ReQ[dQidx(0, l, m)]; + partialImQ[dQidx(blockIdx.x, l, m)] = ImQ[dQidx(0, l, m)]; + } + } + } +} + + +__global__ void centerKernel(Real *rho, Real *bounds, Real *dx, int *n, int n_ghost, Real *partialCenter, Real *partialTotrhosq){ + + __shared__ Real bCenter[3 * CENTERTPB]; + __shared__ Real bTotrhosq[CENTERTPB]; + + for ( int i = threadIdx.x * 3; i < ( threadIdx.x + 1 ) * 3; i++ ) bCenter[i] = 0.; + bTotrhosq[threadIdx.x] = 0.; + + int nreal[3], tid[3]; + for ( int i = 0; i < 3; i++ ) nreal[i] = n[i] - 2 * n_ghost; + int nrealcells = nreal[0] * nreal[1] * nreal[2]; + int tid1d = threadIdx.x + blockIdx.x * blockDim.x; + tid[2] = tid1d / ( nreal[0] * nreal[1] ); + tid[1] = ( tid1d - tid[2] * nreal[0] * nreal[1] ) / nreal[0]; + tid[0] = tid1d - tid[2] * nreal[0] * nreal[1] - tid[1] * nreal[0]; + + Real x[3], rhosq; + int fid; + + while ( tid1d < nrealcells ){ + tid[2] = tid1d / ( nreal[0] * nreal[1] ); + tid[1] = ( tid1d - tid[2] * nreal[0] * nreal[1] ) / nreal[0]; + tid[0] = tid1d - tid[2] * nreal[0] * nreal[1] - tid[1] * nreal[0]; + fid = tidFake(tid[0], tid[1], tid[2], n_ghost, n); + + rhosq = rho[fid] * rho[fid]; + + bTotrhosq[threadIdx.x] += rhosq; + for ( int i = 0; i < 3; i++ ){ + x[i] = bounds[i] + dx[i] * ( tid[i] + 0.5); + bCenter[3 * threadIdx.x + i] += x[i] * rhosq; + } + + tid1d += blockDim.x * gridDim.x; + + } + + __syncthreads(); + + int i = blockDim.x / 2; + while ( i > 0 ){ + if ( threadIdx.x < i ){ + for ( int ii = 0; ii < 3; ii++ ) bCenter[3 * threadIdx.x + ii] += bCenter[3 * ( threadIdx.x + i ) + ii]; + bTotrhosq[threadIdx.x] += bTotrhosq[threadIdx.x + i]; + } + __syncthreads(); + i /= 2; + } + + if ( threadIdx.x == 0 ){ + for ( int i = 0; i < 3; i++) partialCenter[3 * blockIdx.x + i] = bCenter[i]; + partialTotrhosq[blockIdx.x] = bTotrhosq[0]; + } + +} + +//TODO: rmpole should be the distance from the center of the expansion to the nearest boundary cell, not from the center of the domain to the nearest boundary cell +void Grid3D::setMoments(){ + +//Get center of the expansion in the CPU to compare results +/* + int id; + Real totrhosqCPU = 0.; + Real centerCPU[3], x[3]; + + for ( int i = 0; i < 3; i++ ) centerCPU[i] = 0.; + + for ( int k = H.n_ghost; k < H.nz - H.n_ghost; k++ ) { + for ( int j = H.n_ghost; j < H.ny - H.n_ghost; j++ ) { + for ( int i = H.n_ghost; i < H.nx - H.n_ghost; i++ ) { + id = i + j*H.nx + k*H.nx*H.ny; + + // // get the centered cell positions at (i,j,k) + Get_Position(i, j, k, &x[0], &x[1], &x[2]); + + totrhosqCPU += C.density[id] * C.density[id]; + for ( int ii = 0; ii < 3; ii++ ) centerCPU[ii] += x[ii] * C.density[id] * C.density[id]; + + } + } + } + + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, &totrhosqCPU, 1, MPI_CHREAL, MPI_SUM, world); + MPI_Allreduce(MPI_IN_PLACE, centerCPU , 3, MPI_CHREAL, MPI_SUM, world); + #endif + + for ( int i = 0; i < 3; i++ ) centerCPU[i] /= totrhosqCPU; + chprintf("CPU center: %.10e, %.10e, %.10e\n", centerCPU[0], centerCPU[1], centerCPU[2]); +*/ + + Grav.AllocateMemoryBoundaries_GPU(); + Grav.CopyDomainPropertiesToGPU(H.bounds_local, H.n_local_real, H.dxi); + +//Find the center of the expansion according to Couch et al. 2013 + centerKernel<<>>(Grav.Poisson_solver.F.density_d, Grav.dev_bounds, Grav.dev_dx, Grav.dev_n, 0, Grav.dev_partialCenter, Grav.dev_partialTotrhosq); + CudaCheckError(); + + CudaSafeCall( cudaMemcpy(Grav.bufferCenter , Grav.dev_partialCenter , 3 * Grav.centerBlocks * sizeof(Real), cudaMemcpyDeviceToHost) ); + CudaSafeCall( cudaMemcpy(Grav.bufferTotrhosq, Grav.dev_partialTotrhosq, Grav.centerBlocks * sizeof(Real), cudaMemcpyDeviceToHost) ); + + Real totrhosq = 0.; + for ( int i = 0; i < Grav.centerBlocks; i++ ){ + totrhosq += Grav.bufferTotrhosq[i]; + } + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, &totrhosq, 1, MPI_CHREAL, MPI_SUM, world); + #endif + + for ( int i = 0; i < 3; i++ ) Grav.center[i] = 0.; + for ( int i = 0; i < Grav.centerBlocks; i++ ){ + for ( int ii = 0; ii < 3; ii++ ) Grav.center[ii] += Grav.bufferCenter[3 * i + ii]; + } + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, Grav.center, 3, MPI_CHREAL, MPI_SUM, world); + #endif + + for ( int i = 0; i < 3; i++ ) Grav.center[i] /= totrhosq; + + #ifdef POISSON_TEST + for ( int i = 0; i < 3; i++ ) Grav.center[i] = 0.; + #endif + + if ( H.n_step > 0) chprintf(" "); + chprintf("Multipole center: %.10e, %.10e, %.10e\n", Grav.center[0], Grav.center[1], Grav.center[2]); + +//Find the multipole moments + CudaSafeCall( cudaMemcpy( Grav.dev_center, Grav.center, 3*sizeof(Real), cudaMemcpyHostToDevice) ); + + QlmKernel<<>>(Grav.Poisson_solver.F.density_d, Grav.dev_center, Grav.dev_bounds, Grav.dev_dx, H.xdglobal / 2., Grav.dev_n, 0, Grav.dev_partialReQ, Grav.dev_partialImQ); + CudaCheckError(); + + CudaSafeCall( cudaMemcpy(Grav.bufferReQ, Grav.dev_partialReQ, sizeof(Real) * Grav.Qblocks * (1 + LMAX ) * (2 + LMAX ) / 2, cudaMemcpyDeviceToHost) ); + CudaSafeCall( cudaMemcpy(Grav.bufferImQ, Grav.dev_partialImQ, sizeof(Real) * Grav.Qblocks * (1 + LMAX ) * (2 + LMAX ) / 2, cudaMemcpyDeviceToHost) ); + + + for ( int i = 0; i < ( 1 + LMAX ) * ( 2 + LMAX ) / 2; i++ ){ + Grav.ReQ[i] = 0.; + Grav.ImQ[i] = 0.; + } + + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + for ( int b = 0; b < Grav.Qblocks; b++ ){ + + Grav.ReQ[Grav.Qidx(0,l,m)] += Grav.bufferReQ[Grav.Qidx(b, l, m)]; + Grav.ImQ[Grav.Qidx(0,l,m)] += Grav.bufferImQ[Grav.Qidx(b, l, m)]; + + } + + Grav.ReQ[Grav.Qidx(0,l,m)] *= H.dV; + Grav.ImQ[Grav.Qidx(0,l,m)] *= H.dV; + + } + } + + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, Grav.ReQ, (1 + LMAX ) * (2 + LMAX ) / 2, MPI_CHREAL, MPI_SUM, world); + MPI_Allreduce(MPI_IN_PLACE, Grav.ImQ, (1 + LMAX ) * (2 + LMAX ) / 2, MPI_CHREAL, MPI_SUM, world); + #endif + + #ifdef POISSON_TEST + int lmidx; + for ( int l = 0; l <= LMAX; l++ ){ + for ( int m = 0; m <= l; m++ ){ + lmidx = Grav.Qidx(0,l,m); + + chprintf("ReQ[%i][%i]=%.20e\n", l, m, Grav.ReQ[lmidx]); + chprintf("ImQ[%i][%i]=%.20e\n", l, m, Grav.ImQ[lmidx]); + + } + } + #endif + + Grav.FreeMemoryBoundaries_GPU(); + +} + +#endif diff --git a/src/gravity/potential_SOR_3D.cpp b/src/gravity/potential_SOR_3D.cpp index b525388ba..06c80e934 100644 --- a/src/gravity/potential_SOR_3D.cpp +++ b/src/gravity/potential_SOR_3D.cpp @@ -49,10 +49,16 @@ void Potential_SOR_3D::Initialize( Real Lx, Real Ly, Real Lz, Real x_min, Real y TRANSFER_POISSON_BOUNDARIES = false; - chprintf( " Using Poisson Solver: SOR\n"); - chprintf( " SOR: L[ %f %f %f ] N[ %d %d %d ] dx[ %f %f %f ]\n", Lbox_x, Lbox_y, Lbox_z, nx_local, ny_local, nz_local, dx, dy, dz ); + chprintf( " Poisson solver: SOR\n"); + chprintf( " Convergence epsilon: %.5e\n", SOREPSILON); - chprintf( " SOR: Allocating memory...\n"); + #if defined POISSON_TEST || defined TIDES + chprintf( " Maximum angular order: %i\n", LMAX); + #endif + +// chprintf( " SOR: L[ %f %f %f ] N[ %d %d %d ] dx[ %f %f %f ]\n", Lbox_x, Lbox_y, Lbox_z, nx_local, ny_local, nz_local, dx, dy, dz ); + + chprintf( " Allocating memory...\n"); AllocateMemory_CPU(); AllocateMemory_GPU(); @@ -70,7 +76,7 @@ void Potential_SOR_3D::AllocateMemory_CPU( void ){ void Potential_SOR_3D::AllocateMemory_GPU( void ){ - Allocate_Array_GPU_Real( &F.input_d, n_cells_local ); +// Allocate_Array_GPU_Real( &F.input_d, n_cells_local ); Allocate_Array_GPU_Real( &F.density_d, n_cells_local ); Allocate_Array_GPU_Real( &F.potential_d, n_cells_potential ); Allocate_Array_GPU_bool( &F.converged_d, 1 ); @@ -104,7 +110,7 @@ void Potential_SOR_3D::AllocateMemory_GPU( void ){ #endif } - +/* void Potential_SOR_3D::Copy_Input_And_Initialize( Real *input_density, Real Grav_Constant, Real dens_avrg, Real current_a ){ Copy_Input( n_cells_local, F.input_d, input_density, Grav_Constant, dens_avrg, current_a ); @@ -114,7 +120,7 @@ void Potential_SOR_3D::Copy_Input_And_Initialize( Real *input_density, Real Grav potential_initialized = true; } } - +*/ void Potential_SOR_3D::Poisson_Partial_Iteration( int n_step, Real omega, Real epsilon ){ if (n_step == 0 ) Poisson_iteration_Patial_1( n_cells_local, nx_local, ny_local, nz_local, n_ghost, dx, dy, dz, omega, epsilon, F.density_d, F.potential_d, F.converged_h, F.converged_d ); @@ -124,19 +130,18 @@ void Potential_SOR_3D::Poisson_Partial_Iteration( int n_step, Real omega, Real e void Grid3D::Get_Potential_SOR( Real Grav_Constant, Real dens_avrg, Real current_a, struct parameters *P ){ +/* Grav.Poisson_solver.Copy_Input_And_Initialize( Grav.F.density_h, Grav_Constant, dens_avrg, current_a ); + #if defined POISSON_TEST || defined TIDES + setMoments(); + #endif +*/ + //Set Isolated Boundary Conditions Grav.Copy_Isolated_Boundaries_To_GPU( P ); Grav.Poisson_solver.Set_Isolated_Boundary_Conditions( Grav.boundary_flags, P ); - - #if defined POISSON_TEST || defined TIDES - Real epsilon = 1.e-10; - #else - Real epsilon = 1.e-4; - #endif -// chprintf("SOR convergence epsilon: %.5e", epsilon); int max_iter = 10000000; int n_iter = 0; @@ -168,7 +173,7 @@ void Grid3D::Get_Potential_SOR( Real Grav_Constant, Real dens_avrg, Real current Grav.Poisson_solver.TRANSFER_POISSON_BOUNDARIES = false; } - Grav.Poisson_solver.Poisson_Partial_Iteration( 0, omega, epsilon ); + Grav.Poisson_solver.Poisson_Partial_Iteration( 0, omega, SOREPSILON ); if ( set_boundaries ){ Grav.Poisson_solver.TRANSFER_POISSON_BOUNDARIES = true; @@ -176,7 +181,7 @@ void Grid3D::Get_Potential_SOR( Real Grav_Constant, Real dens_avrg, Real current Grav.Poisson_solver.TRANSFER_POISSON_BOUNDARIES = false; } - Grav.Poisson_solver.Poisson_Partial_Iteration( 1, omega, epsilon ); + Grav.Poisson_solver.Poisson_Partial_Iteration( 1, omega, SOREPSILON ); n_iter += 1; @@ -218,7 +223,7 @@ void Potential_SOR_3D::Set_Isolated_Boundary_Conditions( int *boundary_flags, st if ( P->xl_bcnd != 3 && P->xu_bcnd != 3 && P->yl_bcnd != 3 && P->yu_bcnd != 3 && P->zl_bcnd != 3 && P->zu_bcnd != 3 ) return; - chprintf( " Setting Isolated Boundaries \n"); +// chprintf( " Setting Isolated Boundaries \n"); if ( boundary_flags[0] == 3 ) Set_Isolated_Boundary_GPU( 0, 0, F.boundary_isolated_x0_d ); if ( boundary_flags[1] == 3 ) Set_Isolated_Boundary_GPU( 0, 1, F.boundary_isolated_x1_d ); if ( boundary_flags[2] == 3 ) Set_Isolated_Boundary_GPU( 1, 0, F.boundary_isolated_y0_d ); @@ -260,7 +265,7 @@ void Potential_SOR_3D::Copy_Poisson_Boundary_Periodic( int direction, int side ) void Potential_SOR_3D::FreeMemory_GPU( void ){ - Free_Array_GPU_Real( F.input_d ); +// Free_Array_GPU_Real( F.input_d ); Free_Array_GPU_Real( F.density_d ); Free_Array_GPU_Real( F.potential_d ); Free_Array_GPU_Real( F.boundaries_buffer_x0_d ); diff --git a/src/gravity/potential_SOR_3D.h b/src/gravity/potential_SOR_3D.h index 8791267c6..a0050b434 100644 --- a/src/gravity/potential_SOR_3D.h +++ b/src/gravity/potential_SOR_3D.h @@ -53,7 +53,7 @@ class Potential_SOR_3D{ Real *output_h; - Real *input_d; +// Real *input_d; // Real *output_d; Real *density_d; Real *potential_d; @@ -96,7 +96,7 @@ class Potential_SOR_3D{ void AllocateMemory_GPU( void ); void FreeMemory_GPU( void ); void Reset( void ); - void Copy_Input( int n_cells, Real *input_d, Real *input_density_h, Real Grav_Constant, Real dens_avrg, Real current_a ); +// void Copy_Input( int n_cells, Real *input_d, Real *input_density_h, Real Grav_Constant, Real dens_avrg, Real current_a ); void Copy_Output( Real *output_potential ); void Copy_Potential_From_Host( Real *output_potential ); @@ -114,7 +114,7 @@ class Potential_SOR_3D{ void Initialize_Potential( int nx, int ny, int nz, int n_ghost_potential, Real *potential_d, Real *density_d ); - void Copy_Input_And_Initialize( Real *input_density, Real Grav_Constant, Real dens_avrg, Real current_a ); +// void Copy_Input_And_Initialize( Real *input_density, Real Grav_Constant, Real dens_avrg, Real current_a ); void Poisson_iteration( int n_cells, int nx, int ny, int nz, int n_ghost_potential, Real dx, Real dy, Real dz, Real omega, Real epsilon, Real *density_d, Real *potential_d, bool *converged_h, bool *converged_d ); void Poisson_iteration_Patial_1( int n_cells, int nx, int ny, int nz, int n_ghost_potential, Real dx, Real dy, Real dz, Real omega, Real epsilon, Real *density_d, Real *potential_d, bool *converged_h, bool *converged_d ); @@ -154,10 +154,14 @@ class Potential_SOR_3D{ #ifdef MPI_CHOLLA bool Get_Global_Converged( bool converged_local ); #endif + + void Copy_Density_To_GPU(int n_cells, Real *density_h); + void Convert_Density_To_RHS(Real Grav_Constant); + }; #endif //POTENTIAL_SOR_H -#endif //GRAVITY \ No newline at end of file +#endif //GRAVITY diff --git a/src/gravity/potential_SOR_3D_gpu.cu b/src/gravity/potential_SOR_3D_gpu.cu index b7575e283..1147afc1d 100644 --- a/src/gravity/potential_SOR_3D_gpu.cu +++ b/src/gravity/potential_SOR_3D_gpu.cu @@ -28,6 +28,7 @@ void Potential_SOR_3D::Free_Array_GPU_bool( bool *array_dev ){ CudaCheckError(); } +/* __global__ void Copy_Input_Kernel( int n_cells, Real *input_d, Real *density_d, Real Grav_Constant, Real dens_avrg, Real current_a ){ int tid = threadIdx.x + blockIdx.x * blockDim.x; @@ -36,15 +37,14 @@ __global__ void Copy_Input_Kernel( int n_cells, Real *input_d, Real *density_d, #ifdef COSMOLOGY density_d[tid] = 4 * M_PI * Grav_Constant * ( input_d[tid] - dens_avrg ) / current_a; #else - density_d[tid] = 4 * M_PI * Grav_Constant * input_d[tid]; + density_d[tid] = 4 * M_PI * Grav_Constant * density_d[tid]; #endif // if (tid == 0) printf("dens: %f\n", density_d[tid]); } - void Potential_SOR_3D::Copy_Input( int n_cells, Real *input_d, Real *input_density_h, Real Grav_Constant, Real dens_avrg, Real current_a ){ - cudaMemcpy( input_d, input_density_h, n_cells*sizeof(Real), cudaMemcpyHostToDevice ); - +// cudaMemcpy( input_d, input_density_h, n_cells*sizeof(Real), cudaMemcpyHostToDevice ); + cudaMemcpy( F.density_d, input_density_h, n_cells*sizeof(Real), cudaMemcpyHostToDevice ); // set values for GPU kernels int ngrid = (n_cells_local + TPB_SOR - 1) / TPB_SOR; // number of blocks per 1D grid @@ -54,6 +54,36 @@ void Potential_SOR_3D::Copy_Input( int n_cells, Real *input_d, Real *input_densi Copy_Input_Kernel<<>>( n_cells_local, F.input_d, F.density_d, Grav_Constant, dens_avrg, current_a ); } +*/ + +void Potential_SOR_3D::Copy_Density_To_GPU(int n_cells, Real *density_h){ + CudaSafeCall( cudaMemcpy( F.density_d, density_h, n_cells * sizeof(Real), cudaMemcpyHostToDevice) ); +} + +__global__ void Convert_Density_To_RHS_Kernel( int n_cells, Real *density_d, Real Grav_Constant, Real dens_avrg, Real current_a ){ + + int tid = threadIdx.x + blockIdx.x * blockDim.x; + if ( tid >= n_cells ) return; + + #ifdef COSMOLOGY + density_d[tid] = 4 * M_PI * Grav_Constant * ( density_d[tid] - dens_avrg ) / current_a; + #else + density_d[tid] = 4 * M_PI * Grav_Constant * density_d[tid]; + #endif + // if (tid == 0) printf("dens: %f\n", density_d[tid]); +} + + +void Potential_SOR_3D::Convert_Density_To_RHS(Real Grav_Constant){ + // set values for GPU kernels + int ngrid = (n_cells_local + TPB_SOR - 1) / TPB_SOR; + // number of blocks per 1D grid + dim3 dim1dGrid(ngrid, 1, 1); + // number of threads per 1D block + dim3 dim1dBlock(TPB_SOR, 1, 1); + + Convert_Density_To_RHS_Kernel<<>>( n_cells_local, F.density_d, Grav_Constant, 0, 1 ); +} void Grav3D::Copy_Isolated_Boundary_To_GPU_buffer( Real *isolated_boundary_h, Real *isolated_boundary_d, int boundary_size ){ cudaMemcpy( isolated_boundary_d, isolated_boundary_h, boundary_size*sizeof(Real), cudaMemcpyHostToDevice ); @@ -138,22 +168,23 @@ __global__ void Iteration_Step_SOR( int n_cells, Real *density_d, Real *potentia tid_pot = tid_x + tid_y*nx_pot + tid_z*nx_pot*ny_pot; // //Set neighbors ids - int indx_l, indx_r, indx_d, indx_u, indx_b, indx_t, indx_l2, indx_r2, indx_d2, indx_u2, indx_b2, indx_t2; + int indx_l, indx_r, indx_d, indx_u, indx_b, indx_t; +// int indx_l2, indx_r2, indx_d2, indx_u2, indx_b2, indx_t2; indx_l = tid_x-1; //Left - indx_l2 = tid_x-2; //Two to the left +// indx_l2 = tid_x-2; //Two to the left indx_r = tid_x+1; //Right - indx_r2 = tid_x+2; //Two to the right +// indx_r2 = tid_x+2; //Two to the right indx_d = tid_y-1; //Down - indx_d2 = tid_y-2; //Two down +// indx_d2 = tid_y-2; //Two down indx_u = tid_y+1; //Up - indx_u2 = tid_y+2; //Two up +// indx_u2 = tid_y+2; //Two up indx_b = tid_z-1; //Bottom - indx_b2 = tid_z-2; //Two bottom +// indx_b2 = tid_z-2; //Two bottom indx_t = tid_z+1; //Top - indx_t2 = tid_z+2; //Two top +// indx_t2 = tid_z+2; //Two top //Boundary Conditions are loaded to the potential array, the natural indices work! @@ -173,40 +204,41 @@ __global__ void Iteration_Step_SOR( int n_cells, Real *density_d, Real *potentia // indx_b = tid_z == n_ghost ? tid_z+1 : tid_z-1; //Bottom // indx_t = tid_z == nz_pot-n_ghost-1 ? tid_z-1 : tid_z+1; //Top - Real rho, phi_c, phi_l2, phi_l, phi_r, phi_r2, phi_d2, phi_d, phi_u, phi_u2, phi_b2, phi_b, phi_t, phi_t2, phi_new; + Real rho, phi_c, phi_l, phi_r, phi_d, phi_u, phi_b, phi_t, phi_new; +// Real phi_l2, phi_r2, phi_d2, phi_u2, phi_b2, phi_t2; rho = density_d[tid]; phi_c = potential_d[tid_pot]; - phi_l2 = potential_d[ indx_l2 + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; +// phi_l2 = potential_d[ indx_l2 + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; phi_l = potential_d[ indx_l + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; phi_r = potential_d[ indx_r + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; - phi_r2 = potential_d[ indx_r2 + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; - phi_d2 = potential_d[ tid_x + indx_d2 * nx_pot + tid_z * nx_pot * ny_pot ]; +// phi_r2 = potential_d[ indx_r2 + tid_y * nx_pot + tid_z * nx_pot * ny_pot ]; +// phi_d2 = potential_d[ tid_x + indx_d2 * nx_pot + tid_z * nx_pot * ny_pot ]; phi_d = potential_d[ tid_x + indx_d * nx_pot + tid_z * nx_pot * ny_pot ]; phi_u = potential_d[ tid_x + indx_u * nx_pot + tid_z * nx_pot * ny_pot ]; - phi_u2 = potential_d[ tid_x + indx_u2 * nx_pot + tid_z * nx_pot * ny_pot ]; - phi_b2 = potential_d[ tid_x + tid_y * nx_pot + indx_b2 * nx_pot * ny_pot ]; +// phi_u2 = potential_d[ tid_x + indx_u2 * nx_pot + tid_z * nx_pot * ny_pot ]; +// phi_b2 = potential_d[ tid_x + tid_y * nx_pot + indx_b2 * nx_pot * ny_pot ]; phi_b = potential_d[ tid_x + tid_y * nx_pot + indx_b * nx_pot * ny_pot ]; phi_t = potential_d[ tid_x + tid_y * nx_pot + indx_t * nx_pot * ny_pot ]; - phi_t2 = potential_d[ tid_x + tid_y * nx_pot + indx_t2 * nx_pot * ny_pot ]; +// phi_t2 = potential_d[ tid_x + tid_y * nx_pot + indx_t2 * nx_pot * ny_pot ]; /* - if ( tid < 10 ){ - printf("l2: %f\n", phi_l2); - printf("r2: %f\n", phi_r2); - printf("d2: %f\n", phi_d2); - printf("u2: %f\n", phi_u2); - printf("b2: %f\n", phi_b2); - printf("t2: %f\n", phi_t2); - } + if ( tid < 10 ){ + printf("l2: %f\n", phi_l2); + printf("r2: %f\n", phi_r2); + printf("d2: %f\n", phi_d2); + printf("u2: %f\n", phi_u2); + printf("b2: %f\n", phi_b2); + printf("t2: %f\n", phi_t2); + } */ -// 4th order SOR step +// 4th order SOR step /* - phi_new = (1. - omega) *phi_c + phi_new = (1. - omega) *phi_c + ( omega / 90. ) * ( - phi_l2 + 16. * phi_l + 16. * phi_r - phi_r2 - - phi_d2 + 16. * phi_d + 16. * phi_u - phi_u2 - - phi_b2 + 16. * phi_b + 16. * phi_t - phi_t2 - - 12. * dx * dx * rho - ); + - phi_d2 + 16. * phi_d + 16. * phi_u - phi_u2 + - phi_b2 + 16. * phi_b + 16. * phi_t - phi_t2 + - 12. * dx * dx * rho + ); */ phi_new = (1-omega)*phi_c + omega/6*( phi_l + phi_r + phi_d + phi_u + phi_b + phi_t - dx*dx*rho ); diff --git a/src/grid3D.cpp b/src/grid3D.cpp index c3f2fdccf..96e7a4ed7 100644 --- a/src/grid3D.cpp +++ b/src/grid3D.cpp @@ -64,7 +64,7 @@ Grid3D::Grid3D(void) #ifdef PPMC H.n_ghost=4; #endif //PPMC - + #ifdef GRAVITY H.n_ghost_potential_offset = H.n_ghost - N_GHOST_POTENTIAL; #endif @@ -117,8 +117,8 @@ void Grid3D::Initialize(struct parameters *P) int nz_in = P->nz; // Set the CFL coefficient (a global variable) - C_cfl = 0.3; - C_cfl /= 5.; +//TEMPORARY ON: Lower CFL + C_cfl = 0.05; #ifndef MPI_CHOLLA @@ -227,6 +227,12 @@ void Grid3D::Initialize(struct parameters *P) H.density_floor = 0.0; #endif + #ifdef PRESSURE_FLOOR + H.pressure_floor = PRES_FLOOR; + #else + H.pressure_floor = 0.; + #endif + #ifdef TEMPERATURE_FLOOR H.temperature_floor = TEMP_FLOOR; #else @@ -240,7 +246,6 @@ void Grid3D::Initialize(struct parameters *P) H.Output_Initial = true; - } @@ -275,7 +280,11 @@ void Grid3D::AllocateMemory(void) #ifdef POISSON_TEST C.analyticalPotential = (Real *) malloc(H.n_cells * sizeof(Real)); #endif - + + #ifdef TIDES_OUTPUT_POTENTIAL_BH + C.Grav_potential_BH = ( Real *) malloc(H.n_cells * sizeof(Real)); + #endif + #else C.Grav_potential = NULL; #endif @@ -531,9 +540,11 @@ Real Grid3D::Update_Grid(void) // Set the lower limit for density and temperature (Internal Energy) Real U_floor, density_floor; density_floor = H.density_floor; - // Minimum of internal energy from minumum of temperature - U_floor = H.temperature_floor / (gama - 1) / MP * KB * 1e-10;; +// Minimum of internal energy from minumum of temperature. +// To get the minimum U, use the estimated max rho. + U_floor = H.pressure_floor / ( gama - 1. ) / 10.; #ifdef COSMOLOGY + U_floor = H.temperature_floor / (gama - 1) / MP * KB * 1e-10;; U_floor /= Cosmo.v_0_gas * Cosmo.v_0_gas / Cosmo.current_a / Cosmo.current_a; #endif @@ -764,7 +775,15 @@ void Grid3D::FreeMemory(void) #ifdef GRAVITY Grav.FreeMemory_CPU(); #endif - + + #ifdef TIDES + Grav.FreeMemoryBoundaries_GPU(); + #endif + + #ifdef TIDES_OUTPUT_POTENTIAL_BH + free(C.Grav_potential_BH); + #endif + #ifdef PARTICLES Particles.Reset(); #endif diff --git a/src/grid3D.h b/src/grid3D.h index 802568539..da655df3a 100644 --- a/src/grid3D.h +++ b/src/grid3D.h @@ -137,6 +137,8 @@ struct Header * \brief Number of real cells in the z-dimension */ int nz_real; + int n_local_real[3]; + /*! \var xbound */ /* \brief Global domain x-direction minimum */ Real xbound; @@ -185,6 +187,10 @@ struct Header /* \brief Global domain length in z-direction */ Real zdglobal; + /* ! \var blocal */ + /* \brief Local domain minimum values for all coordinates */ + Real bounds_local[3]; + /*! \var dx * \brief x-width of cells */ Real dx; @@ -196,7 +202,11 @@ struct Header /*! \var dz * \brief z-width of cells */ Real dz; - + + Real dxi[3]; + + Real dV; + /*! \var t * \brief Simulation time */ Real t; @@ -224,6 +234,7 @@ struct Header // Values for lower limit for density and temperature Real density_floor; Real temperature_floor; + Real pressure_floor; Real Ekin_avrg; @@ -300,6 +311,7 @@ class Grid3D #ifdef TIDES Star S; #endif + #ifdef COOLING_GRACKLE // Object that contains data for Grackle cooling Cool_GK Cool; @@ -368,6 +380,9 @@ class Grid3D Real *analyticalPotential; #endif + #ifdef TIDES_OUTPUT_POTENTIAL_BH + Real *Grav_potential_BH; + #endif } C; @@ -607,12 +622,14 @@ class Grid3D void Uniform_Grid(); void Zeldovich_Pancake( struct parameters P ); - + + #ifdef TIDES // Initial Conditions for a Polytropic Star void Polytropic_Star( struct parameters &P ); // Relax the polytrope to achive hydrostatic equilibrium void Polytropic_Star_Relaxation( struct parameters &P ); + #endif #ifdef MPI_CHOLLA @@ -746,15 +763,25 @@ class Grid3D #endif #ifdef TIDES - void AccBh(Real posBhx, Real posBhy, Real posBhz, Real *accBhx, Real *accBhy, Real *accBhz); +// void AccBh(Real posBhx, Real posBhy, Real posBhz, Real *accBhx, Real *accBhy, Real *accBhz); void damp(); void updateCOM(); #endif #ifdef POISSON_TEST void poissonTest( struct parameters P ); + int poissonErrorNorm(); #endif + #if defined POISSON_TEST || defined TIDES + void setMoments(); + void setCenter(); + #endif + + #ifdef TIDES_OUTPUT_POTENTIAL_BH + void updatePotBH(); + #endif + }; diff --git a/src/hydro_cuda.cu b/src/hydro_cuda.cu index 644915f4f..021503e78 100644 --- a/src/hydro_cuda.cu +++ b/src/hydro_cuda.cu @@ -565,9 +565,18 @@ __global__ void Calc_dt_3D(Real *dev_conserved, int nx, int ny, int nz, int n_gh max_dti[tid] = fmax((fabs(vx)+cs)/dx, (fabs(vy)+cs)/dy); max_dti[tid] = fmax(max_dti[tid], (fabs(vz)+cs)/dz); max_dti[tid] = fmax(max_dti[tid], 0.0); + + } + #endif + + #ifdef TIDES +// If density is very low (background), basically ignore the Courant condition for the cell +// TODO: Change so that there's not a fixed threshold value + if ( d < 1.e-10 ){ + max_dti[tid] /= 10; } - #endif + } __syncthreads(); diff --git a/src/initial_conditions.cpp b/src/initial_conditions.cpp index 2837e81ba..b9813a313 100644 --- a/src/initial_conditions.cpp +++ b/src/initial_conditions.cpp @@ -31,6 +31,10 @@ void Grid3D::Set_Initial_Conditions(parameters P) { Set_Domain_Properties(P); Set_Gammas(P.gamma); + #ifdef TIDES + S.initialize(P, H.t, H.dt, H.nx, H.ny, H.nz); + #endif + if (strcmp(P.init, "Constant")==0) { Constant(P.rho, P.vx, P.vy, P.vz, P.P); } else if (strcmp(P.init, "Sound_Wave")==0) { @@ -77,19 +81,16 @@ void Grid3D::Set_Initial_Conditions(parameters P) { } else if (strcmp(P.init, "Zeldovich_Pancake")==0) { Zeldovich_Pancake(P); } - #ifdef TIDES else if (strcmp(P.init, "Polytropic_Star")==0) { Polytropic_Star(P); } - #endif//TIDES - + #endif #ifdef POISSON_TEST else if (strcmp(P.init, "poissonTest") == 0) { poissonTest(P); } #endif - else { chprintf ("ABORT: %s: Unknown initial conditions!\n", P.init); chexit(-1); @@ -147,6 +148,19 @@ void Grid3D::Set_Domain_Properties(struct parameters P) H.ydglobal = H.domlen_y; H.zdglobal = H.domlen_z; + H.dV = H.dx * H.dy * H.dz; + H.bounds_local[0] = H.xbound; + H.bounds_local[1] = H.ybound; + H.bounds_local[2] = H.zbound; + + H.dxi[0] = H.dx; + H.dxi[1] = H.dy; + H.dxi[2] = H.dz; + + H.n_local_real[0] = H.nx - 2 * H.n_ghost; + H.n_local_real[1] = H.ny - 2 * H.n_ghost; + H.n_local_real[2] = H.nz - 2 * H.n_ghost; + #else /*MPI_CHOLLA*/ /* set the local domains on each process */ @@ -1305,18 +1319,18 @@ void Grid3D::poissonTest( struct parameters P ){ // Roseanne's density field if ( r < 1. ){ for (int l = 0; l < 6; l++){ - C.density[id] += P.c[l] * pow(r, l) * pow(1. - r * r, 3.) * gsl_sf_legendre_Pl(l, cos(coords[P.d[l]] / r)); + C.density[id] += P.c[l] * pow(r, l) * pow(1. - r * r, 3.) * gsl_sf_legendre_Pl(l, coords[P.d[l]] / r); C.analyticalPotential[id] += P.c[l] * M_PI * ( - 0.5 * pow(r, l + 8.) / ( 2. * l + 9. ) + 2. * pow(r, l + 6.) / ( 2. * l + 7. ) - 3. * pow(r, l + 4.) / ( 2. * l + 5. ) + 2. * pow(r, l + 2.) / ( 2. * l + 3. ) - 0.5 * pow(r, l ) / ( 2. * l + 1. ) - ) * gsl_sf_legendre_Pl(l, cos(coords[P.d[l]] / r)); + ) * gsl_sf_legendre_Pl(l, coords[P.d[l]] / r); } } else{ for ( int l = 0; l < 6; l++){ - C.analyticalPotential[id] += - P.c[l] * 64. * M_PI * gsl_sf_legendre_Pl(l, cos(coords[P.d[l]] / r)) * 3. / ( 2 * l + 9 ) / ( 2 * l + 7 ) / ( 2 * l + 5 ) / ( 2 * l + 3 ) / ( 2 * l + 1 ) / pow(r, l + 1.); + C.analyticalPotential[id] += - P.c[l] * 64. * M_PI * gsl_sf_legendre_Pl(l, coords[P.d[l]] / r) * 3. / ( 2 * l + 9 ) / ( 2 * l + 7 ) / ( 2 * l + 5 ) / ( 2 * l + 3 ) / ( 2 * l + 1 ) / pow(r, l + 1.); } } diff --git a/src/io.cpp b/src/io.cpp index 7d63703b3..fe7ede5e1 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -36,12 +36,12 @@ void Create_Log_File( struct parameters P ){ #endif string file_name ( LOG_FILE_NAME ); - chprintf( "\nCreating Log File: %s \n\n", file_name.c_str() ); + chprintf( "\nCreating Log File: %s\n", file_name.c_str() ); bool file_exists = false; if (FILE *file = fopen(file_name.c_str(), "r")){ file_exists = true; - chprintf( " File exists, appending values: %s \n\n", file_name.c_str() ); + chprintf( " File exists, appending values: %s \n", file_name.c_str() ); fclose( file ); } @@ -55,17 +55,37 @@ void Create_Log_File( struct parameters P ){ out_file << "\n"; out_file << "Run date: " << dt; out_file.close(); - + +//If we're doing tides, create another file where we put the coordinates of the COM and BH + #if defined TIDES && defined OUTPUT_ALWAYS_COM + file_name = ( "orbit_evolution.log" ); + chprintf("Creating Log File: %s\n", file_name.c_str() ); + + file_exists = false; + if (FILE *file = fopen(file_name.c_str(), "r")){ + file_exists = true; + chprintf( " File exists, appending values: %s. Remember to clean repeated lines.\n\n", file_name.c_str() ); + fclose( file ); + } + else{ + out_file.open(file_name.c_str(), ios::app); +// Spaces required so they're centered with %17.10e+space = 18 chars. 5+8+5 for star, 6+6+6 for BH, 4+9+5 for frame + out_file << " t xstar[0] xstar[1] xstar[2] vstar[0] vstar[1] vstar[2] xBH[0] xBH[1] xBH[2] vBH[0] vBH[1] vBH[2] xFrame[0] xFrame[1] xFrame[2] vFrame[0] vFrame[1] vFrame[2] aFrame[0] aFrame[1] aFrame[2]\n"; + out_file.close(); + } + #endif + + chprintf("\n"); + } -void Write_Message_To_Log_File( const char* message ){ +void Write_Message_To_Log_File(string file_name, const char* message ){ #ifdef MPI_CHOLLA if ( procID != 0 ) return; #endif - - string file_name ( LOG_FILE_NAME ); +// string file_name ( LOG_FILE_NAME ); ofstream out_file; out_file.open(file_name.c_str(), ios::app); out_file << message << endl; @@ -452,17 +472,22 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) attribute_id = H5Acreate(file_id, "n_fields", H5T_STD_I32BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_INT, &H.n_fields); status = H5Aclose(attribute_id); - + +//Save some other useful info + attribute_id = H5Acreate(file_id, "CFL", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &C_cfl); + status = H5Aclose(attribute_id); + #ifdef TIDES attribute_id = H5Acreate(file_id, "Mstar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.Mstar); status = H5Aclose(attribute_id); - + attribute_id = H5Acreate(file_id, "Rstar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.Rstar); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "polyN", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "npoly", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.polyN); status = H5Aclose(attribute_id); @@ -473,7 +498,7 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) attribute_id = H5Acreate(file_id, "eta0", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.eta0); status = H5Aclose(attribute_id); - + attribute_id = H5Acreate(file_id, "eta", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.eta); status = H5Aclose(attribute_id); @@ -490,7 +515,7 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.rp); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "tdynStar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "tdynstar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.tdynStar); status = H5Aclose(attribute_id); @@ -501,7 +526,7 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) attribute_id = H5Acreate(file_id, "Mbox", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &S.Mbox); status = H5Aclose(attribute_id); - #endif//TIDES + #endif #ifdef COSMOLOGY attribute_id = H5Acreate(file_id, "H0", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); @@ -594,44 +619,39 @@ void Grid3D::Write_Header_HDF5(hid_t file_id) status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, Real_data); status = H5Aclose(attribute_id); - #ifdef TIDES - - attribute_id = H5Acreate(file_id, "posFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + #ifdef TIDES + attribute_id = H5Acreate(file_id, "xFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.posFrame); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "velFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "vFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.velFrame); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "accFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "aFrame", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.accFrame); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "posBh", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "xBH", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.posBh); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "velBh", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "vBH", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.velBh); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "accBh", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + attribute_id = H5Acreate(file_id, "aBH", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.accBh); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "posSt", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); - status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.posSt); + attribute_id = H5Acreate(file_id, "xstar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.xstar); status = H5Aclose(attribute_id); - attribute_id = H5Acreate(file_id, "velSt", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); - status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.velSt); + attribute_id = H5Acreate(file_id, "vstar", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.vstar); status = H5Aclose(attribute_id); - -// attribute_id = H5Acreate(file_id, "accSt", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); -// status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, S.accSt); -// status = H5Aclose(attribute_id); - #endif//TIDES + #endif // Close the dataspace status = H5Sclose(dataspace_id); @@ -1230,26 +1250,6 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id) dims[2] = nz_dset; dataspace_id = H5Screate_simple(3, dims, NULL); - #ifdef POISSON_TEST - // Copy the analytical potential array to the memory buffer - for (int k=H.n_ghost; k - +#include /* Write the data */ void WriteData(Grid3D &G, struct parameters P, int nfile); @@ -26,6 +26,6 @@ int chprintf(const char * __restrict sdata, ...); void Create_Log_File( struct parameters P ); -void Write_Message_To_Log_File( const char* message ); +void Write_Message_To_Log_File(std::string file_name, const char* message ) ; #endif /*IO_CHOLLA_H*/ diff --git a/src/main.cpp b/src/main.cpp index 85e1c7854..224ccdec7 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -14,6 +14,10 @@ #include "io.h" #include "error_handling.h" +#ifdef PROFILING +#include +#endif + #define OUTPUT //#define CPU_TIME @@ -37,6 +41,8 @@ int main(int argc, char *argv[]) InitializeChollaMPI(&argc, &argv); #endif /*MPI_CHOLLA*/ + printMemoryUsageGPU(); + Real dti = 0; // inverse time step, 1.0 / dt // input parameter variables @@ -45,7 +51,6 @@ int main(int argc, char *argv[]) int nfile = 0; // number of output files Real outtime = 0; // current output time - // read in command line arguments if (argc != 2) { @@ -58,23 +63,22 @@ int main(int argc, char *argv[]) // create the grid Grid3D G; + // read in the parameters parse_params (param_file, &P); // and output to screen - chprintf ("Parameter values: nx = %d, ny = %d, nz = %d, tout = %f, init = %s, boundaries = %d %d %d %d %d %d\n", - P.nx, P.ny, P.nz, P.tout, P.init, P.xl_bcnd, P.xu_bcnd, P.yl_bcnd, P.yu_bcnd, P.zl_bcnd, P.zu_bcnd); - if (strcmp(P.init, "Read_Grid") == 0 ) chprintf ("Input directory: %s\n", P.indir); - chprintf ("Output directory: %s\n", P.outdir); + chprintf ("Parameter values:\n n: [%d, %d, %d]\n Boundaries: %i %i %i %i %i %i\n Gas gamma: %.5e\n Initial conditions: %s\n Final time: %.5e\n", P.nx, P.ny, P.nz, P.xl_bcnd, P.xu_bcnd, P.yl_bcnd, P.yu_bcnd, P.zl_bcnd, P.zu_bcnd, P.gamma, P.init, P.tout); + if (strcmp(P.init, "Read_Grid") == 0 ) chprintf (" Input directory: %s\n", P.indir); + chprintf (" Output directory: %s\n", P.outdir); //Create a Log file to output run-time messages Create_Log_File(P); // initialize the grid G.Initialize(&P); - chprintf("Local number of grid cells: %d %d %d %d\n", G.H.nx_real, G.H.ny_real, G.H.nz_real, G.H.n_cells); // Set initial conditions and calculate first dt - chprintf("Setting initial conditions...\n"); + chprintf("\nSetting initial conditions...\n"); G.Set_Initial_Conditions(P); chprintf("Initial conditions set.\n"); // set main variables for Read_Grid inital conditions @@ -84,12 +88,10 @@ int main(int argc, char *argv[]) nfile = P.nfile*P.nfull; } - #ifdef TIDES - G.S.initialize(P, G.H.t, G.H.dt); -// chprintf("About to update COM\n"); - G.updateCOM(); - #endif - + printHydroParams(); +// chprintf("Local number of grid cells: %d %d %d %d\n", G.H.nx_real, G.H.ny_real, G.H.nz_real, G.H.n_cells); +// chprintf("Local dx: %.5e, %.5e, %.5e\n", G.H.dx, G.H.dy, G.H.dz); + #ifdef DE chprintf("\nUsing Dual Energy Formalism:\n eta_1: %0.3f eta_2: %0.4f\n", DE_ETA_1, DE_ETA_2 ); char *message = (char*)malloc(50 * sizeof(char)); @@ -97,7 +99,6 @@ int main(int argc, char *argv[]) Write_Message_To_Log_File( message ); #endif - #ifdef CPU_TIME G.Timer.Initialize(); #endif @@ -123,8 +124,14 @@ int main(int argc, char *argv[]) G.Compute_Gravitational_Potential( &P); #endif + #ifdef POISSON_TEST + int passed = G.poissonErrorNorm(); + if (passed) exit(0); + else exit(1); + #endif + // Set boundary conditions (assign appropriate values to ghost cells) for hydro and potential - chprintf("Setting boundary conditions...\n"); + chprintf("\nSetting boundary conditions...\n"); G.Set_Boundary_Conditions_Grid(P); chprintf("Boundary conditions set.\n"); @@ -133,44 +140,36 @@ int main(int argc, char *argv[]) G.Get_Particles_Acceleration(); #endif - chprintf("Dimensions of each cell: dx = %f dy = %f dz = %f\n", G.H.dx, G.H.dy, G.H.dz); - chprintf("Ratio of specific heats gamma = %f\n",gama); - chprintf("Nstep = %d Timestep = %f Simulation time = %f\n", G.H.n_step, G.H.dt, G.H.t); + chprintf("\ndx: [%.5e, %.5e, %.5e]\n", G.H.dx, G.H.dy, G.H.dz); + chprintf("\nNstep = %d Timestep = %f Simulation time = %f\n", G.H.n_step, G.H.dt, G.H.t); -/* #ifdef TIDES - //If solving a polytropic star, do the relaxation step to achive hydrostactic equilibrium - P.nfile = nfile; - if (strcmp(P.init, "Polytropic_Star") == 0 && G.S.tRelax > 0.) G.Polytropic_Star_Relaxation( P ); - nfile = P.nfile; - chprintf("nfile after relaxation: %i\n", P.nfile); + if ( strcmp(P.init, "Polytropic_Star") == 0 && G.S.tRelax > 0. ){ + P.nfile = nfile; +// If solving a polytropic star, do the relaxation step to achive hydrostactic equilibrium + G.Polytropic_Star_Relaxation( P ); + nfile = P.nfile; + chprintf("nfile after relaxation: %i\n", P.nfile); + } + G.S.relaxed = 1; #endif -*/ #ifdef OUTPUT if (strcmp(P.init, "Read_Grid") != 0 || G.H.Output_Now ) { + #ifdef TIDES_OUTPUT_POTENTIAL_BH + G.updatePotBH(); + #endif + #ifdef TIDES + G.updateCOM(); + #endif // write the initial conditions to file - chprintf("Writing initial conditions to file...\n"); - WriteData(G, P, nfile); + chprintf("\nWriting initial conditions to file...\n"); + WriteData(G, P, nfile); } // add one to the output file count nfile++; #endif //OUTPUT - #ifdef POISSON_TEST - #ifdef MPI_CHOLLA - MPI_Barrier(world); - #endif//MPI_CHOLLA - chprintf("Poisson equation solved. Exiting now...\n"); - G.Reset(); - - #ifdef MPI_CHOLLA - MPI_Finalize(); - #endif//MPI_CHOLLA - - return 0; - #endif//POISSON_TEST - // increment the next output time outtime += P.outstep; @@ -181,14 +180,20 @@ int main(int argc, char *argv[]) init_min = ReduceRealMin(init); init_max = ReduceRealMax(init); init_avg = ReduceRealAvg(init); - chprintf("Init min: %9.4f max: %9.4f avg: %9.4f\n", init_min, init_max, init_avg); + chprintf("\nInit min: %9.4f max: %9.4f avg: %9.4f\n", init_min, init_max, init_avg); #else printf("Init %9.4f\n", init); #endif //MPI_CHOLLA #endif //CPU_TIME // Evolve the grid, one timestep at a time - chprintf("Starting calculations.\n"); + chprintf("\nStarting calculations.\n\n"); + + #ifdef PROFILING + cudaProfilerStart(); + #endif + + while (G.H.t < P.tout) { chprintf("n_step: %d \n", G.H.n_step + 1 ); @@ -206,19 +211,22 @@ int main(int argc, char *argv[]) //Transfer the particles that moved outside the local domain G.Transfer_Particles_Boundaries(P); #endif - - #ifdef TIDES - G.S.update(G.H.t, G.H.dt); - G.updateCOM(); - #endif + + #if defined TIDES + G.S.update(G.H.t, G.H.dt); + #ifdef OUTPUT_ALWAYS_COM + if ( G.H.t > 0) G.updateCOM(); + #endif + #endif // Advance the grid by one timestep dti = G.Update_Hydro_Grid(); - #ifdef TIDES -// Damp very low densities by a constant factor - G.damp(); - #endif//TIDES + #ifdef TIDES +// TEMPORARY ON: No tides damping +// Damp very low densities by a constant factor +// G.damp(); + #endif // update the simulation time ( t += dt ) G.Update_Time(); @@ -263,6 +271,16 @@ int main(int argc, char *argv[]) if (G.H.t == outtime || G.H.Output_Now ) { #ifdef OUTPUT + + #ifdef TIDES_OUTPUT_POTENTIAL_BH + G.updatePotBH(); + #endif + + #ifdef TIDES + #ifndef OUTPUT_ALWAYS_COM + G.updateCOM(); + #endif + #endif /*output the grid data*/ WriteData(G, P, nfile); // add one to the output file count @@ -291,6 +309,9 @@ int main(int argc, char *argv[]) } /*end loop over timesteps*/ + #ifdef PROFILING + cudaProfilerStop(); + #endif #ifdef CPU_TIME // Print timing statistics @@ -298,7 +319,10 @@ int main(int argc, char *argv[]) G.Timer.Print_Average_Times( P ); #endif - Write_Message_To_Log_File( "Run completed successfully!"); + Write_Message_To_Log_File(LOG_FILE_NAME, "Run completed successfully!"); + +//TEMPORARY ON: Exit with exit(0); + exit(0); // free the grid G.Reset(); diff --git a/src/mpi_routines.cpp b/src/mpi_routines.cpp index b08d22dd1..dc13c4c16 100644 --- a/src/mpi_routines.cpp +++ b/src/mpi_routines.cpp @@ -21,11 +21,13 @@ MPI_Comm world; /*global communicator*/ MPI_Comm node; /*global communicator*/ MPI_Datatype MPI_CHREAL; /*set equal to MPI_FLOAT or MPI_DOUBLE*/ +MPI_Datatype MPI_CHCOMPLEX; #ifdef PARTICLES MPI_Datatype MPI_PART_INT; /*set equal to MPI_INT or MPI_LONG*/ #endif + //MPI_Requests for nonblocking comm MPI_Request *send_request; MPI_Request *recv_request; @@ -133,7 +135,7 @@ void InitializeChollaMPI(int *pargc, char **pargv[]) /*set process ids in comm world*/ MPI_Comm_rank(MPI_COMM_WORLD, &procID); - + /*find number of processes in comm world*/ MPI_Comm_size(MPI_COMM_WORLD, &nproc); @@ -153,7 +155,14 @@ void InitializeChollaMPI(int *pargc, char **pargv[]) #if PRECISION == 2 MPI_CHREAL = MPI_DOUBLE; #endif /*PRECISION*/ - + + #if PRECISION == 1 + MPI_CHCOMPLEX = MPI_COMPLEX; + #endif + #if PRECISION == 2 + MPI_CHCOMPLEX = MPI_DOUBLE_COMPLEX; + #endif + #ifdef PARTICLES #ifdef PARTICLES_LONG_INTS MPI_PART_INT = MPI_LONG; @@ -274,6 +283,10 @@ void DomainDecomposition(struct parameters *P, struct Header *H, int nx_gin, int else H->nz = nz_local+2*H->n_ghost; H->nz_real = nz_local; + H->n_local_real[0] = nx_local; + H->n_local_real[1] = ny_local; + H->n_local_real[2] = nz_local; + // set total number of cells H->n_cells = H->nx * H->ny * H->nz; @@ -742,16 +755,30 @@ Real ReduceRealAvg(Real x) Real ReduceRealSum(Real x) { - Real in = x; - Real out; - Real y; + Real in = x; + Real out; + Real y; - MPI_Allreduce(&in, &out, 1, MPI_CHREAL, MPI_SUM, world); - y = (Real) out; + MPI_Allreduce(&in, &out, 1, MPI_CHREAL, MPI_SUM, world); + y = (Real) out; - return y; + return y; } +/* +std::complex ReduceComplexSum(std::complex x){ + std::complex in = x; + std::complex out; + std::complex y; + + MPI_Allreduce(&in, &out, 1, MPI_CHCOMPLEX, MPI_SUM, world); + + y = (std::complex) out; + + return y; +} +*/ + #ifdef PARTICLES /* MPI reduction wrapper for sum(part_int)*/ Real ReducePartIntSum(part_int_t x) @@ -863,6 +890,10 @@ void Set_Parallel_Domain(Real xmin_global, Real ymin_global, Real zmin_global, R H->yblocal = ymin_local; H->zblocal = zmin_local; + H->bounds_local[0] = xmin_local; + H->bounds_local[1] = ymin_local; + H->bounds_local[2] = zmin_local; + //printf("ProcessID: %d xbound: %f xdglobal: %f xblocal: %f\n", procID, H->xbound, H->xdglobal, H->xblocal); /*perform 1-D first*/ @@ -907,6 +938,11 @@ void Set_Parallel_Domain(Real xmin_global, Real ymin_global, Real zmin_global, R H->dz = H->domlen_z / (H->nz - 2*H->n_ghost); } + H->dxi[0] = H->dx; + H->dxi[1] = H->dy; + H->dxi[2] = H->dz; + H->dV = H->dx * H->dy * H->dz; + /* make sure the domain is properly set for this decomposition*/ if(pd_flag==0) { @@ -1219,13 +1255,13 @@ void TileBlockDecomposition(void) //initialize np_x, np_y, np_z int np_x = 1; - int np_y = 1; - int np_z = 1; - //printf("nproc %d n_gpf %d\n",nproc,n_gpf); + int np_y = 1; + int np_z = 1; + //printf("nproc %d n_gpf %d\n",nproc,n_gpf); - /*find the greatest prime factor of the number of MPI processes*/ + /*find the greatest prime factor of the number of MPI processes*/ n_gpf = greatest_prime_factor(nproc); - //printf("nproc %d n_gpf %d\n",nproc,n_gpf); + //printf("nproc %d n_gpf %d\n",nproc,n_gpf); /*base decomposition on whether n_gpf==2*/ if(n_gpf!=2) @@ -1250,17 +1286,17 @@ void TileBlockDecomposition(void) /*increase ny, nz round-robin*/ while(np_x*np_y*np_z < nproc) { - np_y*=2; - if(np_x*np_y*np_z==nproc) - break; - np_z*=2; + np_y*=2; + if(np_x*np_y*np_z==nproc) + break; + np_z*=2; } } } }else{ - /*nproc is a power of 2*/ + /*nproc is a power of 2*/ /*if we are dealing with two dimensions, we can just assign domain*/ if(nz_global==1) { @@ -1290,21 +1326,21 @@ void TileBlockDecomposition(void) int n_tmp; if(np_z>np_y) { - n_tmp = np_y; - np_y = np_z; - np_z = n_tmp; + n_tmp = np_y; + np_y = np_z; + np_z = n_tmp; } if(np_y>np_x) { - n_tmp = np_x; - np_x = np_y; - np_y = n_tmp; + n_tmp = np_x; + np_x = np_y; + np_y = n_tmp; } if(np_z>np_y) { - n_tmp = np_y; - np_y = np_z; - np_z = n_tmp; + n_tmp = np_y; + np_y = np_z; + np_z = n_tmp; } //save result diff --git a/src/mpi_routines.h b/src/mpi_routines.h index c7a7769e7..39bfdf9e0 100644 --- a/src/mpi_routines.h +++ b/src/mpi_routines.h @@ -18,8 +18,8 @@ extern int root; /*rank of root process*/ extern int procID_node; /*process rank on node*/ extern int nproc_node; /*number of MPI processes on node*/ -extern MPI_Comm world; /*global communicator*/ -extern MPI_Comm node; /*communicator for each node*/ +extern MPI_Comm world; /*global communicator*/ +extern MPI_Comm node; /*communicator for each node*/ extern MPI_Datatype MPI_CHREAL; /*data type describing float precision*/ @@ -39,8 +39,8 @@ extern int source[6]; /* Decomposition flag */ extern int flag_decomp; -#define SLAB_DECOMP 1 //slab decomposition flag -#define BLOCK_DECOMP 2 //block decomposition flag +#define SLAB_DECOMP 1 //slab decomposition flag +#define BLOCK_DECOMP 2 //block decomposition flag //Communication buffers // For SLAB @@ -151,6 +151,9 @@ Real ReduceRealAvg(Real x); /* MPI reduction wrapper for sum(Real)*/ Real ReduceRealSum(Real x); +/* MPI reduction wrapper for sum(Complex)*/ +//std::complex ReduceComplexSum(std::complex x); + #ifdef PARTICLES /* MPI reduction wrapper for sum(part_int)*/ Real ReducePartIntSum(part_int_t x); diff --git a/src/poisson_test.cpp b/src/poisson_test.cpp new file mode 100644 index 000000000..19331ac6d --- /dev/null +++ b/src/poisson_test.cpp @@ -0,0 +1,68 @@ +#ifdef POISSON_TEST +#include "grid3D.h" +#include "io.h" +#include "math.h" + +int Grid3D::poissonErrorNorm(){ + + Real l2norm; + Real deltasq = 0.; + + #ifndef MPI_CHOLLA + int nx_global = H.nx_real; + int ny_global = H.ny_real; + int nz_global = H.nz_real; + #endif + + int apotidx, potidx; + for ( int k = 0; k < H.nz_real; k++ ){ + for ( int j = 0; j < H.ny_real; j++ ){ + for ( int i = 0; i < H.nx_real; i++ ){ + apotidx = ( i + H.n_ghost ) + ( j + H.n_ghost ) * H.nx + ( k + H.n_ghost ) * H.nx * H.ny; + potidx = (i+N_GHOST_POTENTIAL) + (j+N_GHOST_POTENTIAL)*(Grav.nx_local+2*N_GHOST_POTENTIAL) + (k+N_GHOST_POTENTIAL)*(Grav.nx_local+2*N_GHOST_POTENTIAL)*(Grav.ny_local+2*N_GHOST_POTENTIAL); + +// printf("apot pot %.10e %.10e\n", Grav.F.potential_h[potidx], C.analyticalPotential[apotidx]); + if ( fabs(Grav.F.potential_h[potidx]) < 1.e-10 || fabs(C.analyticalPotential[apotidx]) < 1.e-10 ) chprintf("Potential is wrong...\n"); + + deltasq += pow( C.analyticalPotential[apotidx] - Grav.F.potential_h[potidx], 2. ); + } + } + } + + + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, &deltasq, 1, MPI_CHREAL, MPI_SUM, world); + #endif + + l2norm = sqrt( deltasq / nx_global / ny_global / nz_global ); + chprintf("L2 norm = %.20e\n", l2norm); + + Real correctl2norm; + if (nx_global == 64 ){ + correctl2norm = 0.00012863286755550373; + } + else if ( nx_global == 128 ){ + correctl2norm = 3.209325613406217e-5; + } + else if ( nx_global == 256 ){ + correctl2norm = 8.004489514884087e-6; + } + else if (nx_global == 512 ){ + correctl2norm = 1.9801531450853054e-6; + } + else{ + chprintf("Unsupported resolution!"); + exit(-1); + } + + if ( fabs(l2norm / correctl2norm - 1.) < 1.e-10 ){ + return 1; + } + else{ + return 0; + } + + +} + +#endif diff --git a/src/tides/orbit.cpp b/src/tides/orbit.cpp deleted file mode 100644 index fd28b0d6f..000000000 --- a/src/tides/orbit.cpp +++ /dev/null @@ -1,308 +0,0 @@ -#ifdef TIDES - -#include "../global.h" -#include "../grid3D.h" -#include "tides.h" -#include "../io.h" -#include - -#ifdef MPI_CHOLLA -#include "../mpi_routines.h" -#endif - -void Grid3D::updateCOM(){ - - Real posx, posy, posz; - Real velxStTemp = 0.; - Real velyStTemp = 0.; - Real velzStTemp = 0.; - Real posxStTemp = 0.; - Real posyStTemp = 0.; - Real poszStTemp = 0.; - Real totrhoTemp = 0.; - - #ifdef MPI_CHOLLA - Real totrho; - #endif - - Real rho; - int i, j, k, id; - - for (k=H.n_ghost; k + +#define COMTPB (1024) + +#ifdef MPI_CHOLLA +#include "../mpi_routines.h" +#endif + +__global__ void comKernel(Real *rho, Real *momentum_x, Real *momentum_y, Real *momentum_z, Real *bounds, Real *dx, int *n, int n_ghost, Real *partialxstar, Real *partialvstar){ + + __shared__ Real xstar[COMTPB * 3]; + __shared__ Real vstar[COMTPB * 3]; + + for ( int i = 3 * threadIdx.x; i < 3 * ( threadIdx.x + 1 ); i++ ){ + xstar[i] = 0.; + vstar[i] = 0.; + } + + int nreal[3]; + for ( int i = 0; i < 3; i++ ) nreal[i] = n[i] - 2 * n_ghost; + int nrealcells = nreal[0] * nreal[1] * nreal[2]; + + int tid = threadIdx.x + blockIdx.x * blockDim.x; + int tid_z = tid / ( nreal[0] * nreal[1] ); + int tid_y = ( tid - tid_z * nreal[0] * nreal[1] ) / nreal[0]; + int tid_x = tid - tid_z * nreal[0] * nreal[1] - tid_y * nreal[0]; + + Real x[3]; + int fakeid; + + while ( tid < nrealcells ){ + + tid_z = tid / ( nreal[0] * nreal[1] ); + tid_y = ( tid - tid_z * nreal[0] * nreal[1] ) / nreal[0]; + tid_x = tid - tid_z * nreal[0] * nreal[1] - tid_y * nreal[0]; + fakeid = ( tid_z + n_ghost ) * n[0] * n[1] + ( tid_y + n_ghost ) * n[0] + ( tid_x + n_ghost ); + + x[0] = bounds[0] + dx[0] * ( tid_x + 0.5); + x[1] = bounds[1] + dx[1] * ( tid_y + 0.5); + x[2] = bounds[2] + dx[2] * ( tid_z + 0.5); + +// Position of the center of mass + for ( int ii = 0; ii < 3; ii++ ) xstar[threadIdx.x * 3 + ii] += x[ii] * rho[fakeid]; + +// Velocity of the center of mass + vstar[threadIdx.x * 3 ] += momentum_x[fakeid]; + vstar[threadIdx.x * 3 + 1] += momentum_y[fakeid]; + vstar[threadIdx.x * 3 + 2] += momentum_z[fakeid]; + + tid += blockDim.x * gridDim.x; + } + + __syncthreads(); + + int i = blockDim.x / 2; + while ( i > 0 ){ + if ( threadIdx.x < i){ + for ( int ii = 0; ii < 3; ii++ ){ + xstar[threadIdx.x * 3 + ii] += xstar[( threadIdx.x + i) * 3 + ii]; + vstar[threadIdx.x * 3 + ii] += vstar[( threadIdx.x + i) * 3 + ii]; + } + } + __syncthreads(); + i /= 2; + } + + if ( threadIdx.x == 0 ){ + for ( int ii = 0; ii < 3; ii++){ + partialxstar[3 * blockIdx.x + ii] = xstar[ii]; + partialvstar[3 * blockIdx.x + ii] = vstar[ii]; + } + } + +} + +__global__ void potBHKernel(Real *bounds, Real *dx, Real *xFrame, Real *xBH, Real Mbh, int *n, int n_ghost, Real *potBH){ + + int nreal[3], tid[3], fid; + for ( int i = 0; i < 3; i++ ) nreal[i] = n[i] - 2 * n_ghost; + int nrealcells = nreal[0] * nreal[1] * nreal[2]; + + int tid1d = threadIdx.x + blockIdx.x * blockDim.x; + tid[2] = tid1d / ( nreal[0] * nreal[1] ); + tid[1] = ( tid1d - tid[2] * nreal[0] * nreal[1] ) / nreal[0]; + tid[0] = tid1d - tid[2] * nreal[0] * nreal[1] - tid[1] * nreal[0]; + + Real x[3]; + + while ( tid1d < nrealcells ){ + + tid[2] = tid1d / ( nreal[0] * nreal[1] ); + tid[1] = ( tid1d - tid[2] * nreal[0] * nreal[1] ) / nreal[0]; + tid[0] = tid1d - tid[2] * nreal[0] * nreal[1] - tid[1] * nreal[0]; + fid = ( tid[2] + n_ghost ) * n[0] * n[1] + ( tid[1] + n_ghost ) * n[0] + ( tid[0] + n_ghost ); + + for ( int i = 0; i < 3; i++ ) x[i] = bounds[i] + dx[i] * ( tid[i] + 0.5) + xFrame[i] - xBH[i]; + potBH[fid] = - G_CGS * Mbh / sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); + + tid1d += blockDim.x * gridDim.x; + } + +} + +#ifdef TIDES_OUTPUT_POTENTIAL_BH +void Grid3D::updatePotBH(){ + + Real dx[3], bounds[3]; + int n[3]; + dx[0] = H.dx; dx[1] = H.dy; dx[2] = H.dz; + bounds[0] = H.xblocal; bounds[1] = H.yblocal; bounds[2] = H.zblocal; + n[0] = H.nx; n[1] = H.ny; n[2] = H.nz; + + Real *dev_potBH, *dev_bounds, *dev_dx, *dev_xFrame, *dev_xBH; + int *dev_n; + + CudaSafeCall( cudaMalloc( (void**)&dev_bounds, 3 * sizeof(Real) ) ); + CudaSafeCall( cudaMalloc( (void**)&dev_n , 3 * sizeof(int ) ) ); + CudaSafeCall( cudaMalloc( (void**)&dev_dx , 3 * sizeof(Real) ) ); + CudaSafeCall( cudaMalloc( (void**)&dev_xFrame, 3 * sizeof(Real) ) ); + CudaSafeCall( cudaMalloc( (void**)&dev_xBH , 3 * sizeof(Real) ) ); + CudaSafeCall( cudaMalloc( (void**)&dev_potBH , n[0] * n[1] * n[2] * sizeof(Real) ) ); + + CudaSafeCall( cudaMemcpy( dev_bounds, bounds , 3*sizeof(Real), cudaMemcpyHostToDevice) ); + CudaSafeCall( cudaMemcpy( dev_n , n , 3*sizeof(int ), cudaMemcpyHostToDevice) ); + CudaSafeCall( cudaMemcpy( dev_dx , dx , 3*sizeof(Real), cudaMemcpyHostToDevice) ); + CudaSafeCall( cudaMemcpy( dev_xFrame, S.posFrame, 3*sizeof(Real), cudaMemcpyHostToDevice) ); + CudaSafeCall( cudaMemcpy( dev_xBH , S.posBh , 3*sizeof(Real), cudaMemcpyHostToDevice) ); + + potBHKernel<<>>(dev_bounds, dev_dx, dev_xFrame, dev_xBH, S.Mbh, dev_n, H.n_ghost, dev_potBH); + CudaCheckError(); + + CudaSafeCall( cudaMemcpy(C.Grav_potential_BH, dev_potBH, n[0] * n[1] * n[2] * sizeof(Real), cudaMemcpyDeviceToHost) ); + + cudaFree(dev_n); + cudaFree(dev_bounds); + cudaFree(dev_dx); + cudaFree(dev_xFrame); + cudaFree(dev_xBH); + cudaFree(dev_potBH); + +} +#endif + +void Grid3D::updateCOM(){ + + S.Mbox = Grav.ReQ[0] * sqrt( 4 * M_PI ); + Real totrho = S.Mbox / H.dx / H.dy / H.dz; + + Real dx[3], bounds[3]; + int n[3]; + dx[0] = H.dx; dx[1] = H.dy; dx[2] = H.dz; + bounds[0] = H.xblocal; bounds[1] = H.yblocal; bounds[2] = H.zblocal; + n[0] = H.nx; n[1] = H.ny; n[2] = H.nz; + Real *dev_rho, *dev_momentum_x, *dev_momentum_y, *dev_momentum_z, *dev_bounds, *dev_dx, *dev_partialxstar, *dev_partialvstar; + int *dev_n; + +//Allocate memory in GPU + cudaMalloc( (void**)&dev_rho , n[0] * n[1] * n[2] *sizeof(Real) ); + cudaMalloc( (void**)&dev_momentum_x , n[0] * n[1] * n[2] *sizeof(Real) ); + cudaMalloc( (void**)&dev_momentum_y , n[0] * n[1] * n[2] *sizeof(Real) ); + cudaMalloc( (void**)&dev_momentum_z , n[0] * n[1] * n[2] *sizeof(Real) ); + cudaMalloc( (void**)&dev_bounds , 3 * sizeof(Real) ); + cudaMalloc( (void**)&dev_n , 3 * sizeof(int) ); + cudaMalloc( (void**)&dev_dx , 3 * sizeof(Real) ); + cudaMalloc( (void**)&dev_partialxstar , S.comBlocks * 3 * sizeof(Real) ); + cudaMalloc( (void**)&dev_partialvstar , S.comBlocks * 3 * sizeof(Real) ); + +//Copy inputs to GPU + cudaMemcpy( dev_rho , C.density , n[0] * n[1] * n[2] * sizeof(Real), cudaMemcpyHostToDevice); + cudaMemcpy( dev_momentum_x, C.momentum_x, n[0] * n[1] * n[2] * sizeof(Real), cudaMemcpyHostToDevice); + cudaMemcpy( dev_momentum_y, C.momentum_y, n[0] * n[1] * n[2] * sizeof(Real), cudaMemcpyHostToDevice); + cudaMemcpy( dev_momentum_z, C.momentum_z, n[0] * n[1] * n[2] * sizeof(Real), cudaMemcpyHostToDevice); + cudaMemcpy( dev_bounds , bounds , 3*sizeof(Real), cudaMemcpyHostToDevice); + cudaMemcpy( dev_n , n , 3*sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy( dev_dx , dx , 3*sizeof(Real), cudaMemcpyHostToDevice); + +//Call Kernel + comKernel<<>>(dev_rho, dev_momentum_x, dev_momentum_y, dev_momentum_z, dev_bounds, dev_dx, dev_n, H.n_ghost, dev_partialxstar, dev_partialvstar); + +//Copy result to CPU + cudaMemcpy(S.bufferxstar, dev_partialxstar, sizeof(Real) * S.comBlocks * 3, cudaMemcpyDeviceToHost); + cudaMemcpy(S.buffervstar, dev_partialvstar, sizeof(Real) * S.comBlocks * 3, cudaMemcpyDeviceToHost); + +//Free GPU + cudaFree(dev_rho); + cudaFree(dev_momentum_x); + cudaFree(dev_momentum_y); + cudaFree(dev_momentum_z); + cudaFree(dev_n); + cudaFree(dev_bounds); + cudaFree(dev_dx); + cudaFree(dev_partialxstar); + cudaFree(dev_partialvstar); + + for ( int ii = 0; ii < 3; ii++ ){ + S.vstar[ii] = 0.; + S.xstar[ii] = 0.; + for ( int i = 0; i < S.comBlocks; i++ ){ + S.xstar[ii] += S.bufferxstar[3*i + ii]; + S.vstar[ii] += S.buffervstar[3*i + ii]; + } + } + + for ( int i = 0; i < 3; i++ ){ + S.xstar[i] /= totrho; + S.vstar[i] /= totrho; + } + + #ifdef MPI_CHOLLA + MPI_Allreduce(MPI_IN_PLACE, S.xstar, 3, MPI_CHREAL, MPI_SUM, world); + MPI_Allreduce(MPI_IN_PLACE, S.vstar, 3, MPI_CHREAL, MPI_SUM, world); + #endif + + #ifdef OUTPUT_ALWAYS_COM +//Write the COM, etc to the logfile + char *message = (char*)malloc(500 * sizeof(char)); +// Column headers: <--t--> <--------xstar--------> <--------vstar--------> <---------xbh---------> <---------vbh---------> <-------xFrame -------> <-------vFrame -------> <-------aFrame -------> + sprintf(message, "%17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e %17.10e", H.t, S.xstar[0], S.xstar[1], S.xstar[2], S.vstar[0], S.vstar[1], S.vstar[2], S.posBh[0], S.posBh[1], S.posBh[2], S.velBh[0], S.velBh[1], S.velBh[2], S.posFrame[0], S.posFrame[1], S.posFrame[2], S.velFrame[0], S.velFrame[1], S.velFrame[2], S.accFrame[0], S.accFrame[1], S.accFrame[2]); + Write_Message_To_Log_File("orbit_evolution.log", message); + free(message); + #endif + +} + +Real Star::geteta(Real t){ + + Real num; + Real den; + Real totdyn = t / tdynOrb; + Real aux = 3. * totdyn + sqrt( 8. + 9. * totdyn * totdyn ); + + num = - 2. + pow( aux , 2./3.); + den = sqrt(2.) * pow(aux, 1./3.); + + return num / den; + +} + +Real Star::getdeta(Real t){ + + Real totdyn = t / tdynOrb; + Real aux = 3. * totdyn + sqrt( 8. + 9. * totdyn * totdyn ); + Real num = 2. + pow(aux, 2./3.); + Real den = sqrt( 16. + 18. * totdyn * totdyn ) * pow(aux, 1./3.); + + return num / den / tdynOrb; + +} + +Real Star::getddeta(Real t){ + + Real totdyn = t / tdynOrb; + Real aux0 = 8. + 9. * totdyn * totdyn; + Real aux1 = 3. * totdyn + sqrt( 8. + 9. * totdyn * totdyn); + + Real prefac = pow(aux1,1./3.); + Real term1 = 2. * ( sqrt(aux0) - 2. * pow(aux1, 1./3.) ); + Real term2 = 3. * totdyn * ( -6. - ( - 3. * totdyn + sqrt(aux0) ) * pow(aux1, 1./3.)); + Real num = prefac * ( term1 + term2 ); + Real den = 2. * sqrt(2.) * pow(aux0, 3./2.); + + return num / den / tdynOrb / tdynOrb; + +} + +// Updates the coordinates of the frame for t and t + dt / 2 +void Star::updateFrameCoords(Real t, Real dt){ + + Real deta = getdeta (t); + Real ddeta = getddeta(t); + Real MstFrac = Mbh / ( Mstar + Mbh ); + + Real etaExt = geteta (t + dt / 2.); + Real detaExt = getdeta (t + dt / 2.); + Real ddetaExt = getddeta(t + dt / 2.); + + posFrame[0] = MstFrac * rp * ( 1. - eta * eta ); + posFrame[1] = MstFrac * rp * 2. * eta; + posFrame[2] = 0.; + + velFrame[0] = MstFrac * (-2.) * rp * eta * deta; + velFrame[1] = MstFrac * 2. * rp * deta; + velFrame[2] = 0.; + + accFrame[0] = MstFrac * (-2.) * rp * ( deta * deta + eta * ddeta ); + accFrame[1] = MstFrac * 2. * rp * ddeta; + accFrame[2] = 0.; + + posFrameExt[0] = MstFrac * rp * ( 1. - etaExt * etaExt ); + posFrameExt[1] = MstFrac * rp * 2. * etaExt; + posFrameExt[2] = 0.; + + velFrameExt[0] = MstFrac * (-2.) * rp * etaExt * detaExt; + velFrameExt[1] = MstFrac * 2. * rp * detaExt; + velFrameExt[2] = 0.; + + accFrameExt[0] = MstFrac * (-2.) * rp * ( detaExt * detaExt + etaExt * ddetaExt ); + accFrameExt[1] = MstFrac * 2. * rp * ddetaExt; + accFrameExt[2] = 0.; + +} + +void Star::updateBhCoords(Real t, Real dt){ + +// Case with hardcoded BH trajectory + + Real deta = getdeta (t); + Real ddeta = getddeta(t); + Real etaExt = geteta (t + dt / 2.); + Real detaExt = getdeta (t + dt / 2.); + Real ddetaExt = getddeta(t + dt / 2.); + Real MbhFrac = Mstar / ( Mstar + Mbh ); + + posBh[0] = - MbhFrac * rp * ( 1. - eta * eta ); + posBh[1] = - MbhFrac * rp * 2. * eta; + posBh[2] = 0.; + + velBh[0] = - MbhFrac * (-2.) * rp * eta * deta; + velBh[1] = - MbhFrac * 2. * rp * deta; + velBh[2] = 0.; + + accBh[0] = - MbhFrac * (-2.) * rp * ( deta * deta + eta * ddeta ) ; + accBh[1] = - MbhFrac * 2. * rp * ddeta; + accBh[2] = 0.; + + posBhExt[0] = - MbhFrac * rp * ( 1. - etaExt * etaExt ); + posBhExt[1] = - MbhFrac * rp * 2. * etaExt; + posBhExt[2] = 0.; + + velBhExt[0] = - MbhFrac * (-2.) * rp * etaExt * detaExt; + velBhExt[1] = - MbhFrac * 2. * rp * detaExt; + velBhExt[2] = 0.; + + accBhExt[0] = - MbhFrac * (-2.) * rp * ( detaExt * detaExt + etaExt * ddetaExt ) ; + accBhExt[1] = - MbhFrac * 2. * rp * ddetaExt; + accBhExt[2] = 0.; + +} + +#endif diff --git a/src/tides/polytrope_functions.cpp b/src/tides/polytrope_functions.cpp index 8da0a1eb8..215b15ea1 100644 --- a/src/tides/polytrope_functions.cpp +++ b/src/tides/polytrope_functions.cpp @@ -70,14 +70,14 @@ int Binary_Search( int N, Real val, Real *data, int indx_l, int indx_r ){ Real Interpolate( int n, int rootIdx, Real xi, Real *xiVals, Real *thetaVals, Real *dthetaVals ){ // if ( x <= xVals[0] ) return thetaVals[0]; // if ( x >= xVals[n-1] ) return thetaVals[n-1]; - if ( xi < 0. ) chprintf("Error: radius requested < 0" ); -// if ( r > P.Rstar ) chprintf("Error: radius requested > Rstar"); - if ( xi > xiVals[rootIdx] ){ - Real val = thetaVals[rootIdx] + ( xi - xiVals[rootIdx] ) * dthetaVals[rootIdx]; - chprintf("This cell is after the last xi. Returning %.10e\n", val); - return val; - } - if ( xi > xiVals[rootIdx+1]) chprintf("wtf\n"); + if ( xi < 0. ) chprintf("Error: radius requested < 0" ); +// if ( r > P.Rstar ) chprintf("Error: radius requested > Rstar"); + if ( xi > xiVals[rootIdx] ){ + Real val = thetaVals[rootIdx] + ( xi - xiVals[rootIdx] ) * dthetaVals[rootIdx]; + chprintf("This cell is after the last xi. Returning %.10e\n", val); + return val; + } + if ( xi > xiVals[rootIdx+1]) chprintf("wtf\n"); // Find the closest index for which xVal is less than x; int indx = Binary_Search( n, xi, xiVals, 0, n-1 ); @@ -93,14 +93,12 @@ Real Interpolate( int n, int rootIdx, Real xi, Real *xiVals, Real *thetaVals, Re } void Grid3D::Polytropic_Star( struct parameters &P ){ - - chprintf(" Polytrope mass: %.5e g\n", P.Mstar); - chprintf(" Polytrope radius: %.5e cm\n", P.Rstar); - chprintf(" Polytropic index: %.1f\n", P.polyN ); - + + chprintf(" Lane-Emden solver:\n"); + //Solve Lane–Emden equation for the polytrope // int n_points = 500000000; - int n_points = 10000000; + int n_points = 10000000; Real *xi_vals = new Real[n_points]; Real *theta_vals = new Real[n_points]; Real *theta_deriv = new Real[n_points]; @@ -111,22 +109,22 @@ void Grid3D::Polytropic_Star( struct parameters &P ){ xi_min = 1.e-10; - if ( P.polyN == 0. ){ - xi_max = 2.46; - } - else if ( P.polyN == 1. ){ - xi_max = 3.15; - } - else if ( P.polyN == 1.5 ){ - xi_max = 3.65376; - } - else{ - xi_max = 7.; - } + if ( P.polyN == 0. ){ + xi_max = 2.46; + } + else if ( P.polyN == 1. ){ + xi_max = 3.15; + } + else if ( P.polyN == 1.5 ){ + xi_max = 3.65376; + } + else{ + xi_max = 7.; + } dxi = xi_max / ( n_points - 1. ); - xi_vals[0] = 0.; - xi_vals[1] = xi_min; + xi_vals[0] = 0.; + xi_vals[1] = xi_min; for ( int i = 2; i < n_points; i++){ xi_vals[i] = i * dxi; } @@ -135,21 +133,21 @@ void Grid3D::Polytropic_Star( struct parameters &P ){ vector poly_coords; // The first elements of this vector will be the known boundary conditions - poly_coords.push_back( Real2(1., 0.) ); - theta_vals[0] = 1.; - theta_deriv[0] = 0.; + poly_coords.push_back( Real2(1., 0.) ); + theta_vals[0] = 1.; + theta_deriv[0] = 0.; // We can't start the integration from the previous point because there'll be division by zero. Instead, integrate the first point from the known Taylor series solution to the equation - Real thetaTaylor, dthetaTaylor; + Real thetaTaylor, dthetaTaylor; thetaTaylor = 1. - (1./6.) * xi_min * xi_min + P.polyN * pow(xi_min, 4.) / 120. - P.polyN * ( 8. * P.polyN - 5.) * pow(xi_min, 6.) / 15120.; dthetaTaylor = - xi_min / 3. + P.polyN * pow(xi_min, 3.) / 30. - pow(xi_min, 5.) * P.polyN * ( -5. + 8 * P.polyN ) / 2520.; - poly_coords.push_back( Real2( thetaTaylor, dthetaTaylor ) ); + poly_coords.push_back( Real2( thetaTaylor, dthetaTaylor ) ); -// chprintf("Taylor series: %.10e, %.10e", thetaTaylor, dthetaTaylor); +// chprintf("Taylor series: %.10e, %.10e", thetaTaylor, dthetaTaylor); - theta_vals[1] = poly_coords[1].x; + theta_vals[1] = poly_coords[1].x; theta_deriv[1] = poly_coords[1].y; //Solve the polytrope equation using the RK4 module @@ -174,20 +172,20 @@ void Grid3D::Polytropic_Star( struct parameters &P ){ } /* - for ( int i = 0; i < root_indx + 1; i++){ - chprintf("xi = %.10e, theta = %.10e, dtheta = %.10e\n", xi_vals[i], theta_vals[i], theta_deriv[i]); - } + for ( int i = 0; i < root_indx + 1; i++){ + chprintf("xi = %.10e, theta = %.10e, dtheta = %.10e\n", xi_vals[i], theta_vals[i], theta_deriv[i]); + } */ -// Linear interpolation estimate of the root - Real xi_root = ( xi_vals[root_indx + 1] * theta_vals[root_indx] - xi_vals[root_indx] * theta_vals[root_indx + 1] ) / ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); - chprintf( " Root at xi = %.5e. Theta values before and after: %.5e %.5e\n", xi_root, theta_vals[root_indx], theta_vals[root_indx+1] ); +// Linear interpolation estimate of the root + Real xi_root = ( xi_vals[root_indx + 1] * theta_vals[root_indx] - xi_vals[root_indx] * theta_vals[root_indx + 1] ) / ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); + chprintf( " Root at xi = %.5e. Theta before and after: %.5e %.5e\n", xi_root, theta_vals[root_indx], theta_vals[root_indx+1] ); -// Linear extrapolation estimate of the derivative evaluated at the root +// Linear extrapolation estimate of the derivative evaluated at the root Real theta_deriv_root = xi_vals[root_indx + 1] * theta_vals[root_indx] * ( theta_deriv[root_indx - 1] - theta_deriv[root_indx] ); - theta_deriv_root += xi_vals[root_indx - 1] * theta_deriv[root_indx] * ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); - theta_deriv_root += xi_vals[root_indx] * ( theta_vals[root_indx + 1] * theta_deriv[root_indx] - theta_vals[root_indx] * theta_deriv[root_indx - 1] ); - theta_deriv_root /= ( xi_vals[root_indx - 1] - xi_vals[root_indx] ) * ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); + theta_deriv_root += xi_vals[root_indx - 1] * theta_deriv[root_indx] * ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); + theta_deriv_root += xi_vals[root_indx] * ( theta_vals[root_indx + 1] * theta_deriv[root_indx] - theta_vals[root_indx] * theta_deriv[root_indx - 1] ); + theta_deriv_root /= ( xi_vals[root_indx - 1] - xi_vals[root_indx] ) * ( theta_vals[root_indx] - theta_vals[root_indx + 1] ); chprintf( " d(theta)/d(xi) at the root: %.5e\n", theta_deriv_root ); @@ -199,10 +197,10 @@ void Grid3D::Polytropic_Star( struct parameters &P ){ Real K = pressure_central * pow( dens_central, -(P.polyN+1)/P.polyN ); Real alpha = sqrt( (P.polyN + 1) * K / ( 4 * M_PI * G_CGS ) ) * pow( dens_central, (1.-P.polyN)/(2*P.polyN) ); - chprintf( " dens_central / dens_avrg: %f\n", dens_central / dens_avrg ); - chprintf( " pressure_central: %.5e erg/cm^3\n", pressure_central ); + chprintf( " rho_c / rho_av: %.5e g/cm^3\n", dens_central / dens_avrg ); + chprintf( " p_c : %.5e erg/cm^3\n", pressure_central ); Real cs_center = sqrt( pressure_central / dens_central * P.gamma ); - chprintf( " Sound crossing time: %.5e s\n ", P.Rstar / cs_center); + chprintf( " t_cross : %.5e s\n", P.Rstar / cs_center); // chprintf( " K: %f \n", K ); // chprintf( " alpha: %f \n", alpha ); for ( int i=0; i tRelax - if ( t < S.tRelax || dens < 1.e1 * 1.e-10 ){ - vx *= relaxRate; - vy *= relaxRate; - vy *= relaxRate; - } +// Guillochon+ 2013 relaxation +// The first criterion will apply to all cells when t < tRelax, and only to low density cells when t > tRelax + if ( t < S.tRelax || dens < 1.e1 * 1.e-10 ){ + vx *= relaxRate; + vy *= relaxRate; + vy *= relaxRate; + } - v2 = vx*vx + vy*vy + vz*vz; - v = sqrt( v2 ); + v2 = vx*vx + vy*vy + vz*vz; +// v = sqrt( v2 ); -// Compute the energy with the updated kinetic energy - E = U + 0.5*dens*v2; +// Compute the energy with the updated kinetic energy + E = U + 0.5*dens*v2; -// Save the updated values +// Save the updated values C.momentum_x[id] = dens*vx; C.momentum_y[id] = dens*vy; C.momentum_z[id] = dens*vz; @@ -436,7 +429,7 @@ void Grid3D::damp(){ #ifdef MPI_CHOLLA max_speed_global = ReduceRealMax( max_speed ); #endif - */ + */ } @@ -454,46 +447,43 @@ void Grid3D::Polytropic_Star_Relaxation( struct parameters &P ){ WriteData(*this, P, P.nfile); P.nfile++; - while (H.t < S.tRelax ){ + while (H.t < S.tRelax ){ + + chprintf(" Relaxation n_step: %d\n", n_step + 0 ); + + S.update(H.t, H.dt); + updateCOM(); + // calculate the timestep + set_dt(dti); + + // Advance the grid by one timestep + dti = Update_Hydro_Grid(); + + // update the simulation time ( t += dt ) + Update_Time(); + + // add one to the timestep count + n_step++; + + #ifdef GRAVITY + //Compute Gravitational potential for next step + Compute_Gravitational_Potential( &P); + #endif + + //Include the damping terms in momentum and energy + damp(); + +// TODO: Change from number of steps to time so that it's consistent with the rest of the code + // Output + if (n_step % int(P.outstep) == 0){ + WriteData(*this, P, P.nfile); + P.nfile++; + } + + // set boundary conditions for next time step + Set_Boundary_Conditions_Grid(P); - chprintf(" Relaxation n_step: %d\n", n_step + 1 ); - - - S.update(H.t, H.dt); - updateCOM(); - // calculate the timestep - set_dt(dti); - - // Advance the grid by one timestep - dti = Update_Hydro_Grid(); - - // update the simulation time ( t += dt ) - Update_Time(); - - // add one to the timestep count - n_step++; - - #ifdef GRAVITY - //Compute Gravitational potential for next step - Compute_Gravitational_Potential( &P); - #endif - - //Include the damping terms in momentum and energy - damp(); - - // Output - if (n_step % int(P.outstep) == 0){ - WriteData(*this, P, P.nfile); - P.nfile++; - } - - // set boundary conditions for next time step - Set_Boundary_Conditions_Grid(P); - - chprintf("n_step: %d sim time: %10.7f sim timestep: %7.4e \n\n", n_step, H.t, H.dt); - - //Exit the iteartions if converged - // if ( converged ) break; + chprintf("n_step: %d sim time: %10.7f sim timestep: %7.4e \n\n", n_step, H.t, H.dt); } diff --git a/src/tides/tides.cpp b/src/tides/tides.cpp index d236f1da3..3af9a9411 100644 --- a/src/tides/tides.cpp +++ b/src/tides/tides.cpp @@ -5,16 +5,14 @@ #include "../io.h" // Kronecker delta -int kronDelta(int i, int j){ +Real kronDelta(int i, int j){ - if ( i == j ) return 1; - else{ - return 0; - } + if ( i == j ) return 1.; + else return 0.; } -void Star::initialize(struct parameters &P, Real t, Real dt){ +void Star::initialize(struct parameters P, Real t, Real dt, int nx, int ny, int nz){ Mstar = P.Mstar; Mbh = P.Mbh; @@ -22,164 +20,230 @@ void Star::initialize(struct parameters &P, Real t, Real dt){ polyN = P.polyN; //We currently assume that simulations are restarted after the relaxation - if (strcmp(P.init, "Read_Grid")==0){ - relaxed = 1; - } - else{ - relaxed = 0; - } + if (strcmp(P.init, "Read_Grid")==0){ + relaxed = 1; + } + else{ + relaxed = 0; + } //Mass ratio - q = Mbh / Mstar; + q = Mbh / Mstar; //Reduced mass - mu = Mbh * Mstar / ( Mbh + Mstar ); + mu = Mbh * Mstar / ( Mbh + Mstar ); //Tidal radius - rt = Rstar * pow(q, 1./3.); + rt = Rstar * pow(q, 1./3.); //Pericenter distance - rp = P.rprt * rt; + rp = P.rprt * rt; //Initial eta - eta0 = - sqrt(P.r0rt / P.rprt - 1.); + eta0 = - sqrt(P.r0rt / P.rprt - 1.); //Initial distance - r0 = P.r0rt * rt; + r0 = P.r0rt * rt; //Dynamical time of the star and of the orbit - tdynStar = sqrt( pow(Rstar, 3.) / G_CGS / Mstar ); - tdynOrb = sqrt( pow(rp , 3.) / G_CGS / ( Mstar + Mbh ) ); + tdynStar = sqrt( pow(Rstar, 3.) / G_CGS / Mstar ); + tdynOrb = sqrt( pow(rp , 3.) / G_CGS / ( Mstar + Mbh ) ); //Initial time - t0 = sqrt(2.) * tdynOrb * eta0 * ( 1. + eta0 * eta0 / 3. ); - chprintf("t0 = %.20e\n", t0); - - chprintf(" Dynamical time of the star : %.5e\n", tdynStar); - chprintf(" Dynamical time of the orbit: %.5e\n", tdynOrb ); - chprintf(" Initial eta : %.5e\n", eta0 ); + t0 = sqrt(2.) * tdynOrb * eta0 * ( 1. + eta0 * eta0 / 3. ); // Total energy of the star - E0star = ( G_CGS * Mstar * Mstar / Rstar ) * ( 3. / ( polyN - 5. ) + 1. / ( 5. - polyN ) / ( P.gamma - 1. ) ); + E0star = ( G_CGS * Mstar * Mstar / Rstar ) * ( 3. / ( polyN - 5. ) + 1. / ( 5. - polyN ) / ( P.gamma - 1. ) ); //Relaxation time - tRelax = P.tRelaxtDyn * tdynStar; + tRelax = P.tRelaxtDyn * tdynStar; + relaxRate0 = P.relaxRate0; + relaxRateBkgnd = P.relaxRateBkgnd; + + update(t, dt); + + comBlocks = ceil ( nx * ny * nz / COMTPB ); + bufferxstar = (Real *) malloc( sizeof(Real) * comBlocks * 3); + buffervstar = (Real *) malloc( sizeof(Real) * comBlocks * 3); + + chprintf(" Tidal setup flags:\n"); + chprintf(" Framework : "); + #ifdef TIDES_RELATIVISTIC + chprintf("Relativistic"); + #else + chprintf("Newtonian"); + #endif + chprintf("\n"); + chprintf(" Output COM coords : "); + #ifdef TIDES_OUTPUT_ALWAYS_COM + chprintf("At every step"); + #else + chprintf("On HDF5 only"); + #endif + chprintf("\n"); + chprintf(" Output BH potential: "); + #ifdef TIDES_OUTPUT_POTENTIAL_BH + chprintf("True"); + #else + chprintf("False"); + #endif + chprintf("\n"); + + chprintf(" Star:\n"); + chprintf(" n_poly: %.10e\n", polyN); + chprintf(" Mass : %.10e g\n", Mstar); + chprintf(" Radius: %.10e cm\n", Rstar); + chprintf(" t_dyn : %.10e s\n", tdynStar); + + chprintf(" Orbit:\n"); + chprintf(" Mass ratio : %.10e\n", q); + chprintf(" t_dyn : %.10e s\n", tdynOrb); + chprintf(" Tidal radius : %.10e cm\n", rt ); + chprintf(" Initial dist : %.10e cm\n", r0 ); + chprintf(" Periapsis dist: %.10e cm\n", rp ); + chprintf(" Periapsis time: %.10e s\n", -t0 ); + + if ( tRelax > 0 ) chprintf(" Relaxation enabled. Initial relax rate: %f. Background relax rate: %.f\n", relaxRate0, relaxRateBkgnd); + +} + +void Star::update(Real t, Real dt){ + +//The time along the orbit is different from the hydro time because we relax the star, and because the time along the orbit is measured with t = 0 at periapsis. When the star is not relaxed, we hold the star at the initial position along the orbit (but we compute no tidal forces!). After it's relaxed, we compute the time along the orbit accounting for the offset from the initial conditions. Remember that t0 is negative. + +//When relaxed == 0, the tidal tensors aren't used anyways. + if ( relaxed == 0 ){ + tOrb = t0; + } + else{ + tOrb = t + t0; + } + +// Important: first do frames, then tidal tensors since they depend on the frames! + eta = geteta(tOrb); + updateFrameCoords (tOrb, dt); + updateBhCoords (tOrb, dt); + #ifdef TIDES_RELATIVISTIC + updateTidalTensors(); + #endif + +} - update(t, dt); -// E0orb = getEorb(); + +#ifdef TIDES_RELATIVISTIC +//Given a position and a set of tidal tensors, return the tidal potential +Real Star::getTidalPotential(Real *x, Real argCij[3][3], Real argCijk[3][3][3], Real argCijkl[3][3][3][3]){ + + Real tidalPot = 0.; + for ( int i = 0; i < 3; i++ ){ + for ( int j = 0; j < 3; j++){ + tidalPot += 0.5 * argCij[i][j] * x[i] * x[j]; + for ( int k = 0; k < 3; k++){ + tidalPot += (1./6.) * argCijk[i][j][k] * x[i] * x[j] * x[k]; + for ( int l = 0; l < 3; l++){ + tidalPot += (1./24.) * argCijkl[i][j][k][l] * x[i] * x[j] * x[k] * x[l]; + } + } + } + } + + return tidalPot; } +#else +Real Star::getTidalPotential(Real *x){ + + Real dxaux[3], framePot, globalPot; + for ( int i = 0; i < 3; i++ ) dxaux[i] = posBhExt[i] - posFrameExt[i]; + + framePot = - G_CGS * Mbh * ( x[0] * dxaux[0] + x[1] * dxaux[1] + x[2] * dxaux[2] ) / pow(dxaux[0] * dxaux[0] + dxaux[1] * dxaux[1] + dxaux[2] * dxaux[2], 1.5); + globalPot = - G_CGS * Mbh / sqrt( pow((x[0] - dxaux[0]), 2.) + pow(x[1] - dxaux[1], 2.) + pow(x[2] - dxaux[2], 2.) ); + return globalPot - framePot; -// Given a position and a set of tidal tensors, return the tidal potential -Real Star::getTidalPotential(Real x, Real y, Real z, Real Cij[3][3], Real Cijk[3][3][3], Real Cijkl[3][3][3][3]){ - -//Coordinates where the potential is requested - Real coords[3]; - coords[0] = x; - coords[1] = y; - coords[2] = z; - - Real tidalPot; - -// Using tidal tensors - tidalPot = 0.; - for ( int i = 0; i < 3; i++ ){ - for ( int j = 0; j < 3; j++){ - tidalPot += 0.5 * Cij[i][j] * coords[i] * coords[j]; - for ( int k = 0; k < 3; k++){ - tidalPot += (1./6.) + Cijk[i][j][k] * coords[i] * coords[j] * coords[k]; - for ( int l = 0; l < 3; l++){ - tidalPot += (1./24.) * Cijkl[i][j][k][l] * coords[i] * coords[j] * coords[k] * coords[l]; - } - } - } - } - -// Using the exact Newtonian potential -// Real r0 = sqrt( ); -// tidalPot = - G_CGS * Mbh / rOrb - - return tidalPot; } +#endif +#ifdef TIDES_RELATIVISTIC // Updates the tidal tensors, which only depend on the position of the center of the frame. // Updates tensors for t and t + dt / 2, since the latter will be used in the extrapolated potential. -void Star::updateTidalTensors(Real t, Real dt){ - -//Coordinates - Real r = sqrt( pow( posFrame[0] - posBh[0], 2. ) + pow( posFrame[1] - posBh[1], 2. ) + pow( posFrame[2] - posBh[2], 2. )); - Real r2 = r * r; - Real r3 = r2 * r; - Real r4 = r3 * r; - Real r5 = r4 * r; - - Real rExt = sqrt( pow( posFrameExt[0] - posBhExt[0], 2. ) + pow( posFrameExt[1] - posBhExt[1], 2. ) + pow( posFrameExt[2] - posBhExt[2], 2. )); - Real r2Ext = rExt * rExt; - Real r3Ext = r2Ext * rExt; - Real r4Ext = r3Ext * rExt; - Real r5Ext = r4Ext * rExt; - - for ( int i = 0; i < 3; i++ ){ - for ( int j = 0; j < 3; j++ ){ - -// Quadrupole tensor at t - Cij[i][j] = kronDelta(i, j) - 3. * posFrame[i] * posFrame[j] / r2; - Cij[i][j] *= G_CGS * Mbh / r3; - -// Quadrupole tensor at t + dt/2 - extCij[i][j] = kronDelta(i, j) - 3. * posFrameExt[i] * posFrameExt[j] / r2Ext; - extCij[i][j] *= G_CGS * Mbh / r3Ext; - - for ( int k = 0; k < 3; k++){ - -// Octupole tensor at t - Cijk[i][j][k] = 15. * posFrame[i] * posFrame[j] * posFrame[k] / r3 - - 3. * ( posFrame[i] * kronDelta(j, k) + posFrame[j] * kronDelta(i, k) + posFrame[k] * kronDelta(i, j) ) / r; - Cijk[i][j][k] *= G_CGS * Mbh / r4; - -// Octupole tensor at t + dt / 2 - extCijk[i][j][k] = 15. * posFrameExt[i] * posFrameExt[j] * posFrameExt[k] / r3Ext - - 3. * ( posFrameExt[i] * kronDelta(j, k) + posFrameExt[j] * kronDelta(i, k) + posFrameExt[k] * kronDelta(i, j) ) / rExt; - extCijk[i][j][k] *= G_CGS * Mbh / r4Ext; - - for ( int l = 0; l < 3; l++){ - -// Hexadecapole tensor - Cijkl[i][j][k][l] = - 105. * posFrame[i] * posFrame[j] * posFrame[k] * posFrame[l] / r4 - + 15. * ( kronDelta(i, l) * posFrame[j] * posFrame[k] - + kronDelta(j, l) * posFrame[i] * posFrame[k] - + kronDelta(k, l) * posFrame[i] * posFrame[j] - + kronDelta(i, j) * posFrame[k] * posFrame[l] - + kronDelta(j, k) * posFrame[i] * posFrame[l] - + kronDelta(i, k) * posFrame[j] * posFrame[l] - ) / r2 - - 3. * ( kronDelta(i, j) * kronDelta(k, l) - + kronDelta(j, k) * kronDelta(i, l) - + kronDelta(i, k) * kronDelta(j, l) - ); - Cijkl[i][j][k][l] *= G_CGS * Mbh / r5; - -// Hexadecapole tensor at t + dt / 2 - extCijkl[i][j][k][l] = - 105. * posFrameExt[i] * posFrameExt[j] * posFrameExt[k] * posFrameExt[l] / r4Ext - + 15. * ( kronDelta(i, l) * posFrameExt[j] * posFrameExt[k] - + kronDelta(j, l) * posFrameExt[i] * posFrameExt[k] - + kronDelta(k, l) * posFrameExt[i] * posFrameExt[j] - + kronDelta(i, j) * posFrameExt[k] * posFrameExt[l] - + kronDelta(j, k) * posFrameExt[i] * posFrameExt[l] - + kronDelta(i, k) * posFrameExt[j] * posFrameExt[l] - ) / r2Ext - - 3. * ( kronDelta(i, j) * kronDelta(k, l) - + kronDelta(j, k) * kronDelta(i, l) - + kronDelta(i, k) * kronDelta(j, l) - ); - extCijkl[i][j][k][l] *= G_CGS * Mbh / r5Ext; - - } - } - } - } - +void Star::updateTidalTensors(){ + + Real r2 = pow( posFrame[0] - posBh[0], 2. ) + pow( posFrame[1] - posBh[1], 2. ) + pow( posFrame[2] - posBh[2], 2. ); + Real r = sqrt(r2); + Real r3 = r2 * r; + Real r4 = r3 * r; + Real r5 = r4 * r; + + Real r2Ext = pow( posFrameExt[0] - posBhExt[0], 2. ) + pow( posFrameExt[1] - posBhExt[1], 2. ) + pow( posFrameExt[2] - posBhExt[2], 2. ); + Real rExt = sqrt( r2Ext ); + Real r3Ext = r2Ext * rExt; + Real r4Ext = r3Ext * rExt; + Real r5Ext = r4Ext * rExt; + + Real bigx[3], bigxExt[3]; + for ( int i = 0; i < 3; i++ ) bigx[i] = posFrame[i] - posBh[i]; + for ( int i = 0; i < 3; i++ ) bigxExt[i] = posFrameExt[i] - posBhExt[i]; + + for ( int i = 0; i < 3; i++ ){ + for ( int j = 0; j < 3; j++ ){ + +// Quadrupole tensor at t + Cij[i][j] = kronDelta(i, j) - 3. * bigx[i] * bigx[j] / r2; + Cij[i][j] *= G_CGS * Mbh / r3; + +// Quadrupole tensor at t + dt/2 + extCij[i][j] = kronDelta(i, j) - 3. * bigxExt[i] * bigxExt[j] / r2Ext; + extCij[i][j] *= G_CGS * Mbh / r3Ext; + + for ( int k = 0; k < 3; k++){ + +// Octupole tensor at t + Cijk[i][j][k] = 15. * bigx[i] * bigx[j] * bigx[k] / r3 + - 3. * ( bigx[i] * kronDelta(j, k) + bigx[j] * kronDelta(i, k) + bigx[k] * kronDelta(i, j) ) / r; + Cijk[i][j][k] *= G_CGS * Mbh / r4; + +// Octupole tensor at t + dt / 2 + extCijk[i][j][k] = 15. * bigxExt[i] * bigxExt[j] * bigxExt[k] / r3Ext + - 3. * ( bigxExt[i] * kronDelta(j, k) + bigxExt[j] * kronDelta(i, k) + bigxExt[k] * kronDelta(i, j) ) / rExt; + extCijk[i][j][k] *= G_CGS * Mbh / r4Ext; + + for ( int l = 0; l < 3; l++){ + +// Hexadecapole tensor + Cijkl[i][j][k][l] = - 105. * bigx[i] * bigx[j] * bigx[k] * bigx[l] / r4 + + 15. * ( kronDelta(i, l) * bigx[j] * bigx[k] + + kronDelta(j, l) * bigx[i] * bigx[k] + + kronDelta(k, l) * bigx[i] * bigx[j] + + kronDelta(i, j) * bigx[k] * bigx[l] + + kronDelta(j, k) * bigx[i] * bigx[l] + + kronDelta(i, k) * bigx[j] * bigx[l] + ) / r2 + - 3. * ( kronDelta(i, j) * kronDelta(k, l) + + kronDelta(j, k) * kronDelta(i, l) + + kronDelta(i, k) * kronDelta(j, l) + ); + Cijkl[i][j][k][l] *= G_CGS * Mbh / r5; + +// Hexadecapole tensor at t + dt / 2 + extCijkl[i][j][k][l] = - 105. * bigxExt[i] * bigxExt[j] * bigxExt[k] * bigxExt[l] / r4Ext + + 15. * ( kronDelta(i, l) * bigxExt[j] * bigxExt[k] + + kronDelta(j, l) * bigxExt[i] * bigxExt[k] + + kronDelta(k, l) * bigxExt[i] * bigxExt[j] + + kronDelta(i, j) * bigxExt[k] * bigxExt[l] + + kronDelta(j, k) * bigxExt[i] * bigxExt[l] + + kronDelta(i, k) * bigxExt[j] * bigxExt[l] + ) / r2Ext + - 3. * ( kronDelta(i, j) * kronDelta(k, l) + + kronDelta(j, k) * kronDelta(i, l) + + kronDelta(i, k) * kronDelta(j, l) + ); + extCijkl[i][j][k][l] *= G_CGS * Mbh / r5Ext; + + } + } + } + } } - +#endif #endif diff --git a/src/tides/tides.h b/src/tides/tides.h index 3c014445d..352f9b80e 100644 --- a/src/tides/tides.h +++ b/src/tides/tides.h @@ -6,99 +6,104 @@ #include #include #include "../global.h" -//#include "../grid3D.h" + +#define COMTPB (1024) class Star { public: - //Star - Real Mstar; - Real Rstar; - Real polyN; - Real tRelax; - Real tdynStar; - int relaxed; + Real Mstar; + Real Rstar; + Real polyN; + Real tRelax; + Real relaxRate0; + Real relaxRateBkgnd; + Real tdynStar; + int relaxed; //Total mass in the box. Used to compute mass loss for close encounters - Real Mbox; + Real Mbox; //Tides - Real rp; - Real r0; - Real tdynOrb; - Real Mbh; - Real q; - Real mu; - Real rt; - Real t0; - Real tOrb; - Real eta0; - Real eta; - Real E0orb; - Real Eorb; - Real E0star; - Real DEoE; + Real rp; + Real r0; + Real tdynOrb; + Real Mbh; + Real q; + Real mu; + Real rt; + Real t0; + Real tOrb; + Real eta0; + Real eta; + Real E0orb; + Real Eorb; + Real E0star; + Real DEoE; //Coordinates of the center of the box - Real posFrame[3]; - Real velFrame[3]; - Real accFrame[3]; + Real posFrame[3]; + Real velFrame[3]; + Real accFrame[3]; //Extrapolated - Real posFrameExt[3]; - Real velFrameExt[3]; - Real accFrameExt[3]; + Real posFrameExt[3]; + Real velFrameExt[3]; + Real accFrameExt[3]; -// Coordinates of the bh - Real posBh[3]; - Real velBh[3]; - Real accBh[3]; +//Coordinates of the bh + Real posBh[3]; + Real velBh[3]; + Real accBh[3]; //Extrapolated - Real posBhExt[3]; - Real velBhExt[3]; - Real accBhExt[3]; - -// Coordinates of the star - Real posSt[3]; - Real velSt[3]; - Real accSt[3]; + Real posBhExt[3]; + Real velBhExt[3]; + Real accBhExt[3]; - Real posStExt[3]; - Real velStExt[3]; - Real accStExt[3]; +// Coordinates of the star + Real xstar[3]; + Real vstar[3]; + Real astar[3]; //Tidal tensors at the current time and at t + dt / 2 + #ifdef TIDES_RELATIVISTIC + Real extCij[3][3]; + Real extCijk[3][3][3]; + Real extCijkl[3][3][3][3]; - Real extCij[3][3]; - Real extCijk[3][3][3]; - Real extCijkl[3][3][3][3]; - - Real Cij[3][3]; - Real Cijk[3][3][3]; - Real Cijkl[3][3][3][3]; + Real Cij[3][3]; + Real Cijk[3][3][3]; + Real Cijkl[3][3][3][3]; + #endif //Functions that change the state of S - void initialize(struct parameters &P, Real t, Real dt); - void update(Real t, Real dt); - void updateFrameCoords(Real t, Real dt); - void updateBhCoords(Real t, Real dt); - void updateTidalTensors(Real t, Real dt); + void initialize(struct parameters P, Real t, Real dt, int nx, int ny, int nz); + void update(Real t, Real dt); + void updateFrameCoords(Real t, Real dt); + void updateBhCoords(Real t, Real dt); + #ifdef TIDES_RELATIVISTIC + void updateTidalTensors(); + #endif //Value of eta (proxy for time) and its first two derivatives with respect to time. These are used to track the coordinates of the center of the frame at all times analytically - Real geteta(Real t); - Real getdeta(Real t); - Real getddeta(Real t); + Real geteta(Real t); + Real getdeta(Real t); + Real getddeta(Real t); //Returns the tidal potential given a set of tidal tensors - Real getTidalPotential(Real x, Real y, Real z, Real Cij[3][3], Real Cijk[3][3][3], Real Cijkl[3][3][3][3]); - -//Returns the position of the center of the frame at time t - void getFrameAndBhPos(Real t, Real *xFrame, Real *yFrame, Real *zFrame, Real *xBh, Real *yBh, Real *zBh); - - Real getEorb(); + #ifdef TIDES_RELATIVISTIC + Real getTidalPotential(Real *x, Real argCij[3][3], Real argCijk[3][3][3], Real argCijkl[3][3][3][3]); + #else + Real getTidalPotential(Real *x); + #endif + +//Used for computing the center of mass position and speed in the GPU + int comBlocks; + Real *bufferxstar; + Real *buffervstar; }; diff --git a/tests/poisson_test/CMakeLists.txt b/tests/poisson_test/CMakeLists.txt new file mode 100644 index 000000000..22c2b3b90 --- /dev/null +++ b/tests/poisson_test/CMakeLists.txt @@ -0,0 +1,18 @@ +project(cholla) + +add_executable(cholla ../../src/main.cpp) + + +enable_testing() + +add_test(poisson64 ../../cholla.sor poissonParameterFiles/poisson64.txt) +set_tests_properties(poisson64 PROPERTIES WILL_FAIL FALSE) + +add_test(poisson128 ../../cholla.sor poissonParameterFiles/poisson128.txt) +set_tests_properties(poisson128 PROPERTIES WILL_FAIL FALSE) + +add_test(poisson256 ../../cholla.sor poissonParameterFiles/poisson256.txt) +set_tests_properties(poisson256 PROPERTIES WILL_FAIL FALSE) + +add_test(poisson512 ../../cholla.sor poissonParameterFiles/poisson512.txt) +set_tests_properties(poisson512 PROPERTIES WILL_FAIL FALSE) diff --git a/tests/poisson_test/poissonParameterFiles/poisson128.txt b/tests/poisson_test/poissonParameterFiles/poisson128.txt new file mode 100644 index 000000000..26c2b886f --- /dev/null +++ b/tests/poisson_test/poissonParameterFiles/poisson128.txt @@ -0,0 +1,62 @@ +# +# Parameter File for the 3D Polytropic Star. +# +###################################### + +# number of grid cells in the x dimension +nx=128 + +# number of grid cells in the y dimension +ny=128 + +# number of grid cells in the z dimension +nz=128 + +# output time +tout=0. + +c0=0.75 +c1=0.5 +c2=0.75 +c3=1. +c4=1. +c5=1. + +#c0=1. +#c1=0. +#c2=0. +#c3=0. +#c4=0. +#c5=0. + +d0=0 +d1=0 +d2=0 +d3=0 +d4=0 +d5=0 + +# how often to output +outstep=100000000000 + +# value of gamma +gamma=1.66666667 + +# name of initial conditions +init=poissonTest + +# domain properties +xmin=-2. +ymin=-2. +zmin=-2. +xlen=4. +ylen=4. +zlen=4. + +# type of boundary conditions +xl_bcnd=3 +xu_bcnd=3 +yl_bcnd=3 +yu_bcnd=3 +zl_bcnd=3 +zu_bcnd=3 diff --git a/tests/poisson_test/poissonParameterFiles/poisson256.txt b/tests/poisson_test/poissonParameterFiles/poisson256.txt new file mode 100644 index 000000000..7c78b1ae1 --- /dev/null +++ b/tests/poisson_test/poissonParameterFiles/poisson256.txt @@ -0,0 +1,62 @@ +# +# Parameter File for the 3D Polytropic Star. +# +###################################### + +# number of grid cells in the x dimension +nx=256 + +# number of grid cells in the y dimension +ny=256 + +# number of grid cells in the z dimension +nz=256 + +# output time +tout=0. + +c0=0.75 +c1=0.5 +c2=0.75 +c3=1. +c4=1. +c5=1. + +#c0=1. +#c1=0. +#c2=0. +#c3=0. +#c4=0. +#c5=0. + +d0=0 +d1=0 +d2=0 +d3=0 +d4=0 +d5=0 + +# how often to output +outstep=100000000000 + +# value of gamma +gamma=1.66666667 + +# name of initial conditions +init=poissonTest + +# domain properties +xmin=-2. +ymin=-2. +zmin=-2. +xlen=4. +ylen=4. +zlen=4. + +# type of boundary conditions +xl_bcnd=3 +xu_bcnd=3 +yl_bcnd=3 +yu_bcnd=3 +zl_bcnd=3 +zu_bcnd=3 diff --git a/tests/poisson_test/poissonParameterFiles/poisson512.txt b/tests/poisson_test/poissonParameterFiles/poisson512.txt new file mode 100644 index 000000000..efa7c10fb --- /dev/null +++ b/tests/poisson_test/poissonParameterFiles/poisson512.txt @@ -0,0 +1,62 @@ +# +# Parameter File for the 3D Polytropic Star. +# +###################################### + +# number of grid cells in the x dimension +nx=512 + +# number of grid cells in the y dimension +ny=512 + +# number of grid cells in the z dimension +nz=512 + +# output time +tout=0. + +c0=0.75 +c1=0.5 +c2=0.75 +c3=1. +c4=1. +c5=1. + +#c0=1. +#c1=0. +#c2=0. +#c3=0. +#c4=0. +#c5=0. + +d0=0 +d1=0 +d2=0 +d3=0 +d4=0 +d5=0 + +# how often to output +outstep=100000000000 + +# value of gamma +gamma=1.66666667 + +# name of initial conditions +init=poissonTest + +# domain properties +xmin=-2. +ymin=-2. +zmin=-2. +xlen=4. +ylen=4. +zlen=4. + +# type of boundary conditions +xl_bcnd=3 +xu_bcnd=3 +yl_bcnd=3 +yu_bcnd=3 +zl_bcnd=3 +zu_bcnd=3 diff --git a/tests/poisson_test/poissonParameterFiles/poisson64.txt b/tests/poisson_test/poissonParameterFiles/poisson64.txt new file mode 100644 index 000000000..95bc0412d --- /dev/null +++ b/tests/poisson_test/poissonParameterFiles/poisson64.txt @@ -0,0 +1,62 @@ +# +# Parameter File for the 3D Polytropic Star. +# +###################################### + +# number of grid cells in the x dimension +nx=64 + +# number of grid cells in the y dimension +ny=64 + +# number of grid cells in the z dimension +nz=64 + +# output time +tout=0. + +c0=0.75 +c1=0.5 +c2=0.75 +c3=1. +c4=1. +c5=1. + +#c0=1. +#c1=0. +#c2=0. +#c3=0. +#c4=0. +#c5=0. + +d0=0 +d1=0 +d2=0 +d3=0 +d4=0 +d5=0 + +# how often to output +outstep=100000000000 + +# value of gamma +gamma=1.66666667 + +# name of initial conditions +init=poissonTest + +# domain properties +xmin=-2. +ymin=-2. +zmin=-2. +xlen=4. +ylen=4. +zlen=4. + +# type of boundary conditions +xl_bcnd=3 +xu_bcnd=3 +yl_bcnd=3 +yu_bcnd=3 +zl_bcnd=3 +zu_bcnd=3 diff --git a/tests/poisson_test/run_output.log b/tests/poisson_test/run_output.log new file mode 100644 index 000000000..b71b35319 --- /dev/null +++ b/tests/poisson_test/run_output.log @@ -0,0 +1,46 @@ + +Run date: Mon Apr 5 13:33:14 2021 + +Run date: Mon Apr 5 13:33:14 2021 + +Run date: Mon Apr 5 13:33:16 2021 + +Run date: Mon Apr 5 13:36:17 2021 + +Run date: Mon Apr 5 13:36:17 2021 + +Run date: Mon Apr 5 13:36:19 2021 + +Run date: Mon Apr 5 13:38:25 2021 + +Run date: Mon Apr 5 13:38:26 2021 + +Run date: Mon Apr 5 13:38:27 2021 + +Run date: Mon Apr 5 13:38:47 2021 + +Run date: Mon Apr 5 13:38:48 2021 + +Run date: Mon Apr 5 13:39:05 2021 + +Run date: Mon Apr 5 13:39:06 2021 + +Run date: Mon Apr 5 13:39:07 2021 + +Run date: Mon Apr 5 13:39:17 2021 + +Run date: Mon Apr 5 13:42:33 2021 + +Run date: Mon Apr 5 13:42:33 2021 + +Run date: Mon Apr 5 13:42:35 2021 + +Run date: Mon Apr 5 13:42:45 2021 + +Run date: Mon Apr 5 14:51:54 2021 + +Run date: Mon Apr 5 14:51:54 2021 + +Run date: Mon Apr 5 14:51:56 2021 + +Run date: Mon Apr 5 14:52:06 2021 diff --git a/tests/poisson_test/test_results b/tests/poisson_test/test_results new file mode 100644 index 000000000..736b807a1 --- /dev/null +++ b/tests/poisson_test/test_results @@ -0,0 +1,15 @@ +[HANDLER_OUTPUT] +Test project /home/rcastroy/src/cholla/tests/poisson_test + + Start 1: poisson64 +1/4 Test #1: poisson64 ........................ Passed 0.40 sec + Start 2: poisson128 +2/4 Test #2: poisson128 ....................... Passed 1.48 sec + Start 3: poisson256 +3/4 Test #3: poisson256 ....................... Passed 9.86 sec + Start 4: poisson512 +4/4 Test #4: poisson512 ....................... Passed 81.90 sec + +100% tests passed, 0 tests failed out of 4 + +Total Test time (real) = 93.65 sec