From 7400fa374ff523077023b957af63f9717da4a205 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 14:22:41 -0700 Subject: [PATCH 01/62] First commit adding properties needed for non-standard cosmologies. We need to track explicitly Omega_R (radiation), w0, and wa for time-variable dark energy models. We also need to write to disk the expansion history, so we can verify that it's correct for unusual cosmologies. We've implemented much of that framework here, but will need to test. --- src/analysis/io_analysis.cpp | 24 +++++++++ src/cosmology/cosmology.cpp | 6 ++- src/cosmology/cosmology.h | 7 +++ src/cosmology/cosmology_functions.cpp | 68 +++++++++++++++++++++++-- src/global/global.cpp | 9 +++- src/global/global.h | 4 ++ src/gravity/gravity_functions.cpp | 3 ++ src/io/io.cpp | 4 +- src/particles/particles_3D.h | 6 ++- src/particles/particles_dynamics.cpp | 5 +- src/particles/particles_dynamics_gpu.cu | 28 ++++++---- src/system_tests/system_tester.cpp | 7 +++ 12 files changed, 147 insertions(+), 24 deletions(-) diff --git a/src/analysis/io_analysis.cpp b/src/analysis/io_analysis.cpp index 962503dea..cd1c0c3ff 100644 --- a/src/analysis/io_analysis.cpp +++ b/src/analysis/io_analysis.cpp @@ -73,6 +73,18 @@ void Grid3D::Write_Skewers_Header_HDF5(hid_t file_id) status = H5Aclose(attribute_id); attribute_id = H5Acreate(file_id, "Omega_b", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_b); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "Omega_K", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_K); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "Omega_R", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_R); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "w0", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.w0); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "wa", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.wa); status = H5Aclose(attribute_id); #endif @@ -525,6 +537,18 @@ void Grid3D::Write_Analysis_Header_HDF5(hid_t file_id) attribute_id = H5Acreate(file_id, "Omega_b", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_b); status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "Omega_K", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_K); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "Omega_R", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.Omega_R); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "w0", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.w0); + status = H5Aclose(attribute_id); + attribute_id = H5Acreate(file_id, "wa", H5T_IEEE_F64BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &Cosmo.wa); + status = H5Aclose(attribute_id); #endif status = H5Sclose(dataspace_id); diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index 6575798e2..de2cf9d9d 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -15,8 +15,12 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part H0 /= 1000; //[km/s / kpc] Omega_M = P->Omega_M; Omega_L = P->Omega_L; - Omega_K = 1 - (Omega_M + Omega_L); + Omega_R = P->Omega_R; + Omega_K = 1 - (Omega_M + Omega_L + Omega_R); Omega_b = P->Omega_b; + w0 = P->w0; + wa = P->wa; + if (strcmp(P->init, "Read_Grid") == 0) { // Read scale factor value from Particles diff --git a/src/cosmology/cosmology.h b/src/cosmology/cosmology.h index 1e7c9bd1c..2a9272bb2 100644 --- a/src/cosmology/cosmology.h +++ b/src/cosmology/cosmology.h @@ -19,6 +19,9 @@ class Cosmology Real Omega_L; Real Omega_K; Real Omega_b; + Real Omega_R; + Real w0; + Real wa; Real cosmo_G; Real cosmo_h; @@ -66,6 +69,10 @@ class Cosmology Real Get_da_from_dt(Real dt); Real Get_dt_from_da(Real da); + + //write expansion history log file + void Create_Expansion_History_File(struct Parameters P); + void Write_Expansion_History_Entry(Real t, Real a); }; #endif diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index f00c7e174..6b21fb8ef 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -5,13 +5,19 @@ #include "../grid/grid_enum.h" #include "../io/io.h" + + + void Grid3D::Initialize_Cosmology(struct Parameters *P) { chprintf("Initializing Cosmology... \n"); Cosmo.Initialize(P, Grav, Particles); + //Create expansion history log file + Create_Log_File(P, Cosmo); + // Change to comoving Cosmological System - Change_Cosmological_Frame_Sytem(true); + Change_Cosmological_Frame_System(true); if (fabs(Cosmo.current_a - Cosmo.next_output) < 1e-5) { H.Output_Now = true; @@ -23,14 +29,16 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Real Cosmology::Get_da_from_dt(Real dt) { Real a2 = current_a * current_a; - Real a_dot = sqrt(Omega_M / current_a + a2 * Omega_L + Omega_K) * H0; + Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); + Real a_dot = sqrt(Omega_R/a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return a_dot * dt; } Real Cosmology::Get_dt_from_da(Real da) { Real a2 = current_a * current_a; - Real a_dot = sqrt(Omega_M / current_a + a2 * Omega_L + Omega_K) * H0; + Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); + Real a_dot = sqrt(Omega_R/a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return da / a_dot; } @@ -38,11 +46,13 @@ Real Cosmology::Get_Hubble_Parameter(Real a) { Real a2 = a * a; Real a3 = a2 * a; - Real factor = (Omega_M / a3 + Omega_K / a2 + Omega_L); + Real a4 = a2 * a2; + Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); + Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L*fac_de); return H0 * sqrt(factor); } -void Grid3D::Change_Cosmological_Frame_Sytem(bool forward) +void Grid3D::Change_Cosmological_Frame_System(bool forward) { if (forward) { chprintf(" Converting to Cosmological Comoving System\n"); @@ -130,4 +140,52 @@ void Grid3D::Change_GAS_Frame_System(bool forward) } } + +// writing the expansion history + +void Create_Expansion_History_File(struct Parameters P) +{ + if (not Is_Root_Proc()) { + return; + } + + std::string file_name(EXPANSION_HISTORY_FILE_NAME); + chprintf("\nCreating Expansion History File: %s \n\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()); + fclose(file); + } + + // current date/time based on current system + time_t now = time(0); + // convert now to string form + char *dt = ctime(&now); + + std::string message = std::string(Omega_M) + " " + std::string(Omega_K); + message += " " + std::string(Omega_R) + " " + std::string(Omega_L) + " " + std::string(w0) + " " + std::string(wa); + + std::ofstream out_file; + out_file.open(file_name.c_str(), std::ios::app); + out_file << "\n"; + out_file << "Run date: " << dt; + out_file.close(); +} + +void Write_Expansion_History_Entry(void) +{ + if (not Is_Root_Proc()) { + return; + } + + std::string message = std::string(Cosmo.t_secs/MYR) + " " + std::string(Cosmo.current_a); + std::string file_name(EXPANSION_HISTORY_FILE_NAME); + std::ofstream out_file; + out_file.open(file_name.c_str(), std::ios::app); + out_file << message.c_str() << std::endl; + out_file.close(); +} + #endif diff --git a/src/global/global.cpp b/src/global/global.cpp index 5cce16bb5..28c3d749b 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -104,7 +104,8 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Init_redshift", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", + "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); @@ -368,6 +369,12 @@ bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameter parms->Omega_L = atof(value); } else if (strcmp(name, "Omega_b") == 0) { parms->Omega_b = atof(value); + } else if (strcmp(name, "Omega_R") == 0) { + parms->Omega_R = atof(value); + } else if (strcmp(name, "w0") == 0) { + parms->w0 = atof(value); + } else if (strcmp(name, "wa") == 0) { + parms->wa = atof(value); #endif // COSMOLOGY #ifdef TILED_INITIAL_CONDITIONS } else if (strcmp(name, "tile_length") == 0) { diff --git a/src/global/global.h b/src/global/global.h index 3d83aa246..cea176a95 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -321,10 +321,14 @@ struct Parameters { Real Omega_M; Real Omega_L; Real Omega_b; + Real Omega_R; + Real w0; + Real wa; Real Init_redshift; Real End_redshift; char scale_outputs_file[MAXLEN]; // File for the scale_factor output values // for cosmological simulations +#define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" #endif // COSMOLOGY #ifdef TILED_INITIAL_CONDITIONS Real tile_length; diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 07c2697ea..2c24cc196 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -127,6 +127,9 @@ void Grid3D::set_dt_Gravity() chprintf(" t_physical: %f Myr dt_physical: %f Myr\n", Cosmo.t_secs / MYR, Cosmo.dt_secs / MYR); Particles.dt = dt_physical; + //Write expansion history + Write_Expansion_History_Entry(); + #else // Not Cosmology // If NOT using COSMOLOGY diff --git a/src/io/io.cpp b/src/io/io.cpp index 9e1fbe534..12f18fe9c 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -120,7 +120,7 @@ void Write_Data(Grid3D &G, struct Parameters P, int nfile) #endif #ifdef COSMOLOGY - G.Change_Cosmological_Frame_Sytem(false); + G.Change_Cosmological_Frame_System(false); #endif #ifndef ONLY_PARTICLES @@ -175,7 +175,7 @@ void Write_Data(Grid3D &G, struct Parameters P, int nfile) } else { chprintf(" Saved Snapshot: %d z:%f\n", nfile, G.Cosmo.current_z); } - G.Change_Cosmological_Frame_Sytem(true); + G.Change_Cosmological_Frame_System(true); chprintf("\n"); G.H.Output_Now = false; #endif diff --git a/src/particles/particles_3D.h b/src/particles/particles_3D.h index 5fca6054d..2ddd146e7 100644 --- a/src/particles/particles_3D.h +++ b/src/particles/particles_3D.h @@ -274,13 +274,15 @@ class Particles3D Real *pos_y_dev, Real *pos_z_dev, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, - Real cosmo_h, Real Omega_M, Real Omega_L, Real Omega_K); + Real cosmo_h, Real Omega_M, Real Omega_L, Real Omega_K, + Real Omega_R, Real w0, Real wa); void Advance_Particles_KDK_Step2_GPU_function(part_int_t n_local, Real dt, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev); void Advance_Particles_KDK_Step2_Cosmo_GPU_function(part_int_t n_local, Real delta_a, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, - Real cosmo_h, Real Omega_M, Real Omega_L, Real Omega_K); + Real cosmo_h, Real Omega_M, Real Omega_L, Real Omega_K, + Real Omega_R, Real w0, Real wa); part_int_t Compute_Particles_GPU_Array_Size(part_int_t n); int Select_Particles_to_Transfer_GPU(int direction, int side); void Copy_Transfer_Particles_to_Buffer_GPU(int n_transfer, int direction, int side, Real *send_buffer, diff --git a/src/particles/particles_dynamics.cpp b/src/particles/particles_dynamics.cpp index 39aeba6c7..43e482023 100644 --- a/src/particles/particles_dynamics.cpp +++ b/src/particles/particles_dynamics.cpp @@ -96,7 +96,8 @@ void Grid3D::Advance_Particles_KDK_Step1_GPU() Particles.Advance_Particles_KDK_Step1_Cosmo_GPU_function( Particles.n_local, Cosmo.delta_a, Particles.pos_x_dev, Particles.pos_y_dev, Particles.pos_z_dev, Particles.vel_x_dev, Particles.vel_y_dev, Particles.vel_z_dev, Particles.grav_x_dev, Particles.grav_y_dev, - Particles.grav_z_dev, Cosmo.current_a, Cosmo.H0, Cosmo.cosmo_h, Cosmo.Omega_M, Cosmo.Omega_L, Cosmo.Omega_K); + Particles.grav_z_dev, Cosmo.current_a, Cosmo.H0, Cosmo.cosmo_h, Cosmo.Omega_M, Cosmo.Omega_L, Cosmo.Omega_K, + Cosmo.Omega_R, Cosmo.w0, Cosmo.wa); #else Particles.Advance_Particles_KDK_Step1_GPU_function(Particles.n_local, Particles.dt, Particles.pos_x_dev, Particles.pos_y_dev, Particles.pos_z_dev, Particles.vel_x_dev, @@ -112,7 +113,7 @@ void Grid3D::Advance_Particles_KDK_Step2_GPU() Particles.Advance_Particles_KDK_Step2_Cosmo_GPU_function( Particles.n_local, Cosmo.delta_a, Particles.vel_x_dev, Particles.vel_y_dev, Particles.vel_z_dev, Particles.grav_x_dev, Particles.grav_y_dev, Particles.grav_z_dev, Cosmo.current_a, Cosmo.H0, Cosmo.cosmo_h, - Cosmo.Omega_M, Cosmo.Omega_L, Cosmo.Omega_K); + Cosmo.Omega_M, Cosmo.Omega_L, Cosmo.Omega_K, Cosmo.Omega_R, Cosmo.w0, Cosmo.wa); #else Particles.Advance_Particles_KDK_Step2_GPU_function(Particles.n_local, Particles.dt, Particles.vel_x_dev, Particles.vel_y_dev, Particles.vel_z_dev, Particles.grav_x_dev, diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index 817040dca..a788ee2ff 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -18,11 +18,14 @@ // FUTURE FIX: The Hubble function was defined here because I couldn't get it // form other file, tried -dc flag when compiling buu paris broke. -__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K) +__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, + Real Omega_R, Real w0, Real wa) { Real a2 = a * a; Real a3 = a2 * a; - Real factor = (Omega_M / a3 + Omega_K / a2 + Omega_L); + Real a4 = a2 * a2; + Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); + Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L*fac_de); return H0 * sqrt(factor); } #endif @@ -189,7 +192,8 @@ __global__ void Advance_Particles_KDK_Step1_Cosmo_Kernel(part_int_t n_local, Rea Real *pos_z_dev, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, Real cosmo_h, - Real Omega_M, Real Omega_L, Real Omega_K) + Real Omega_M, Real Omega_L, Real Omega_K, Real Omega_R, + Real w0, Real wa) { part_int_t tid = blockIdx.x * blockDim.x + threadIdx.x; if (tid >= n_local) { @@ -205,8 +209,8 @@ __global__ void Advance_Particles_KDK_Step1_Cosmo_Kernel(part_int_t n_local, Rea da_half = da / 2; a_half = current_a + da_half; - H = Get_Hubble_Parameter_dev(current_a, H0, Omega_M, Omega_L, Omega_K); - H_half = Get_Hubble_Parameter_dev(a_half, H0, Omega_M, Omega_L, Omega_K); + H = Get_Hubble_Parameter_dev(current_a, H0, Omega_M, Omega_L, Omega_K, Omega_R, w0, wa); + H_half = Get_Hubble_Parameter_dev(a_half, H0, Omega_M, Omega_L, Omega_K, Omega_R, w0, wa); dt = da / (current_a * H) * cosmo_h; dt_half = da / (a_half * H_half) * cosmo_h / (a_half); @@ -233,7 +237,8 @@ __global__ void Advance_Particles_KDK_Step1_Cosmo_Kernel(part_int_t n_local, Rea __global__ void Advance_Particles_KDK_Step2_Cosmo_Kernel(part_int_t n_local, Real da, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, Real cosmo_h, - Real Omega_M, Real Omega_L, Real Omega_K) + Real Omega_M, Real Omega_L, Real Omega_K, Real, Omega_R, + Real w0, Real wa) { part_int_t tid = blockIdx.x * blockDim.x + threadIdx.x; if (tid >= n_local) { @@ -249,7 +254,8 @@ __global__ void Advance_Particles_KDK_Step2_Cosmo_Kernel(part_int_t n_local, Rea da_half = da / 2; a_half = current_a - da_half; - dt = da / (current_a * Get_Hubble_Parameter_dev(current_a, H0, Omega_M, Omega_L, Omega_K)) * cosmo_h; + dt = da / (current_a * Get_Hubble_Parameter_dev(current_a, H0, Omega_M, Omega_L, Omega_K, Omega_R, w0, wa)); + dt *= cosmo_h; // Advance velocities by the second half a step vel_x_dev[tid] = (a_half * vel_x + 0.5 * dt * grav_x_dev[tid]) / current_a; @@ -262,7 +268,7 @@ void Particles3D::Advance_Particles_KDK_Step1_Cosmo_GPU_function(part_int_t n_lo Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, Real cosmo_h, Real Omega_M, Real Omega_L, - Real Omega_K) + Real Omega_K, Real Omega_R, Real w0, Real wa) { // set values for GPU kernels int ngrid = (n_local - 1) / TPB_PARTICLES + 1; @@ -275,7 +281,7 @@ void Particles3D::Advance_Particles_KDK_Step1_Cosmo_GPU_function(part_int_t n_lo if (n_local > 0) { hipLaunchKernelGGL(Advance_Particles_KDK_Step1_Cosmo_Kernel, dim1dGrid, dim1dBlock, 0, 0, n_local, delta_a, pos_x_dev, pos_y_dev, pos_z_dev, vel_x_dev, vel_y_dev, vel_z_dev, grav_x_dev, grav_y_dev, - grav_z_dev, current_a, H0, cosmo_h, Omega_M, Omega_L, Omega_K); + grav_z_dev, current_a, H0, cosmo_h, Omega_M, Omega_L, Omega_K, Omega_R, w0, wa); GPU_Error_Check(cudaDeviceSynchronize()); // GPU_Error_Check(); } @@ -285,7 +291,7 @@ void Particles3D::Advance_Particles_KDK_Step2_Cosmo_GPU_function(part_int_t n_lo Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, Real cosmo_h, Real Omega_M, Real Omega_L, - Real Omega_K) + Real Omega_K, Real Omega_R, Real w0, Real wa) { // set values for GPU kernels int ngrid = (n_local - 1) / TPB_PARTICLES + 1; @@ -298,7 +304,7 @@ void Particles3D::Advance_Particles_KDK_Step2_Cosmo_GPU_function(part_int_t n_lo if (n_local > 0) { hipLaunchKernelGGL(Advance_Particles_KDK_Step2_Cosmo_Kernel, dim1dGrid, dim1dBlock, 0, 0, n_local, delta_a, vel_x_dev, vel_y_dev, vel_z_dev, grav_x_dev, grav_y_dev, grav_z_dev, current_a, H0, cosmo_h, - Omega_M, Omega_L, Omega_K); + Omega_M, Omega_L, Omega_K, Omega_R, w0, wa); GPU_Error_Check(cudaDeviceSynchronize()); // GPU_Error_Check(); } diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 1888fd752..9dec001c6 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -353,6 +353,13 @@ void system_test::SystemTestRunner::launchCholla() } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } +#ifdef COSMOLOGY + try { + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); + } catch (const std::filesystem::filesystem_error &error) { + // This file might not exist and isn't required so don't worry if it doesn't exist + } +#endif //COSMOLOGY } // ============================================================================= From edb86eb063e6e990550d78d58732d9c8ed924558 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 15:37:33 -0700 Subject: [PATCH 02/62] Desperate attempt to get clang-tidy to allow my code to be committed. --- src/global/global.cpp | 6 +++--- src/gravity/gravity_functions.cpp | 2 +- src/particles/particles_dynamics_gpu.cu | 8 ++++---- src/system_tests/system_tester.cpp | 5 +++-- 4 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 28c3d749b..01acc5542 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,9 +103,9 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", - "w0", "wa", "Init_redshift", + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", + "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 2c24cc196..bf5f43228 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -20,7 +20,7 @@ #include "../model/disk_galaxy.h" // #endif -// Set delta_t when usi#ng gravity +// Set delta_t when using gravity void Grid3D::set_dt_Gravity() { // Delta_t for the hydro diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index a788ee2ff..fe0caa70b 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -18,14 +18,14 @@ // FUTURE FIX: The Hubble function was defined here because I couldn't get it // form other file, tried -dc flag when compiling buu paris broke. -__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, - Real Omega_R, Real w0, Real wa) +__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, + Real Omega_K, Real Omega_R, Real w0, Real wa) { Real a2 = a * a; Real a3 = a2 * a; Real a4 = a2 * a2; - Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); - Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L*fac_de); + Real fac_de = pow(a, -3 * (1 + w0 + aw)) * exp(-3 * wa * (1-a)); + Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L * fac_de); return H0 * sqrt(factor); } #endif diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 9dec001c6..4e258413a 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -355,11 +355,12 @@ void system_test::SystemTestRunner::launchCholla() } #ifdef COSMOLOGY try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", + _outputDirectory + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } -#endif //COSMOLOGY +#endif } // ============================================================================= From a263dce29ea2d8d27898cca9de14acf455ac5c50 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 15:46:34 -0700 Subject: [PATCH 03/62] More futile attempts at sating clang-tidy --- src/cosmology/cosmology_functions.cpp | 16 ++++++++-------- src/particles/particles_dynamics_gpu.cu | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 6b21fb8ef..cd7bf75ca 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -28,17 +28,17 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Real Cosmology::Get_da_from_dt(Real dt) { - Real a2 = current_a * current_a; - Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); - Real a_dot = sqrt(Omega_R/a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; + Real a2 = current_a * current_a; + Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return a_dot * dt; } Real Cosmology::Get_dt_from_da(Real da) { - Real a2 = current_a * current_a; - Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); - Real a_dot = sqrt(Omega_R/a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; + Real a2 = current_a * current_a; + Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return da / a_dot; } @@ -47,8 +47,8 @@ Real Cosmology::Get_Hubble_Parameter(Real a) Real a2 = a * a; Real a3 = a2 * a; Real a4 = a2 * a2; - Real fac_de = pow(a,-3*(1+w0+aw))*exp(-3*wa*(1-a)); - Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L*fac_de); + Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L * fac_de); return H0 * sqrt(factor); } diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index fe0caa70b..fcbcdf335 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -18,13 +18,13 @@ // FUTURE FIX: The Hubble function was defined here because I couldn't get it // form other file, tried -dc flag when compiling buu paris broke. -__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, +__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, Real Omega_R, Real w0, Real wa) { Real a2 = a * a; Real a3 = a2 * a; Real a4 = a2 * a2; - Real fac_de = pow(a, -3 * (1 + w0 + aw)) * exp(-3 * wa * (1-a)); + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L * fac_de); return H0 * sqrt(factor); } From c7b9ef038b6fc5ce720e47151ec53f5956620c51 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 15:48:41 -0700 Subject: [PATCH 04/62] Clang tidiness issues --- src/cosmology/cosmology_functions.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index cd7bf75ca..432db05df 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -29,7 +29,7 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Real Cosmology::Get_da_from_dt(Real dt) { Real a2 = current_a * current_a; - Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return a_dot * dt; } @@ -37,7 +37,7 @@ Real Cosmology::Get_da_from_dt(Real dt) Real Cosmology::Get_dt_from_da(Real da) { Real a2 = current_a * current_a; - Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return da / a_dot; } @@ -47,7 +47,7 @@ Real Cosmology::Get_Hubble_Parameter(Real a) Real a2 = a * a; Real a3 = a2 * a; Real a4 = a2 * a2; - Real fac_de = pow(a,-3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real factor = (Omega_R / a4 + Omega_M / a3 + Omega_K / a2 + Omega_L * fac_de); return H0 * sqrt(factor); } @@ -180,7 +180,7 @@ void Write_Expansion_History_Entry(void) return; } - std::string message = std::string(Cosmo.t_secs/MYR) + " " + std::string(Cosmo.current_a); + std::string message = std::string(Cosmo.t_secs / MYR) + " " + std::string(Cosmo.current_a); std::string file_name(EXPANSION_HISTORY_FILE_NAME); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); From a5aa1c9db0f73ad1d4647b31720a623da564ecf9 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 15:52:44 -0700 Subject: [PATCH 05/62] Literally who knows. --- src/cosmology/cosmology_functions.cpp | 10 +++++----- src/global/global.h | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 432db05df..dc8100ac4 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -13,7 +13,7 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) chprintf("Initializing Cosmology... \n"); Cosmo.Initialize(P, Grav, Particles); - //Create expansion history log file + // Create expansion history log file Create_Log_File(P, Cosmo); // Change to comoving Cosmological System @@ -28,17 +28,17 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Real Cosmology::Get_da_from_dt(Real dt) { - Real a2 = current_a * current_a; + Real a2 = current_a * current_a; Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); - Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; + Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return a_dot * dt; } Real Cosmology::Get_dt_from_da(Real da) { - Real a2 = current_a * current_a; + Real a2 = current_a * current_a; Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); - Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; + Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return da / a_dot; } diff --git a/src/global/global.h b/src/global/global.h index cea176a95..86cd70460 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -329,7 +329,7 @@ struct Parameters { char scale_outputs_file[MAXLEN]; // File for the scale_factor output values // for cosmological simulations #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" -#endif // COSMOLOGY +#endif /*COSMOLOGY*/ #ifdef TILED_INITIAL_CONDITIONS Real tile_length; #endif // TILED_INITIAL_CONDITIONS From c9fc6d0c52ac70ca4d1609fe858f8b0b969b97d3 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 15:56:38 -0700 Subject: [PATCH 06/62] No idea how to address clang-tidy --- src/cosmology/cosmology.h | 2 +- src/global/global.h | 3 ++- src/gravity/gravity_functions.cpp | 2 +- src/system_tests/system_tester.cpp | 3 +-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cosmology/cosmology.h b/src/cosmology/cosmology.h index 2a9272bb2..53252ab16 100644 --- a/src/cosmology/cosmology.h +++ b/src/cosmology/cosmology.h @@ -70,7 +70,7 @@ class Cosmology Real Get_da_from_dt(Real dt); Real Get_dt_from_da(Real da); - //write expansion history log file + // write expansion history log file void Create_Expansion_History_File(struct Parameters P); void Write_Expansion_History_Entry(Real t, Real a); }; diff --git a/src/global/global.h b/src/global/global.h index 86cd70460..3d73f43a3 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -328,8 +328,9 @@ struct Parameters { Real End_redshift; char scale_outputs_file[MAXLEN]; // File for the scale_factor output values // for cosmological simulations + #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" -#endif /*COSMOLOGY*/ +#endif // COSMOLOGY #ifdef TILED_INITIAL_CONDITIONS Real tile_length; #endif // TILED_INITIAL_CONDITIONS diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index bf5f43228..aea3acb15 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -127,7 +127,7 @@ void Grid3D::set_dt_Gravity() chprintf(" t_physical: %f Myr dt_physical: %f Myr\n", Cosmo.t_secs / MYR, Cosmo.dt_secs / MYR); Particles.dt = dt_physical; - //Write expansion history + // Write expansion history Write_Expansion_History_Entry(); #else // Not Cosmology diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 4e258413a..0abebdeb8 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -355,8 +355,7 @@ void system_test::SystemTestRunner::launchCholla() } #ifdef COSMOLOGY try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", - _outputDirectory + "/expansion_history.txt"); + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } From 6417f3ce2f6e475d6fe11cd4b5cd2bb104074ea1 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:02:42 -0700 Subject: [PATCH 07/62] Even more random space editing for clang-tidy --- src/cosmology/cosmology.cpp | 1 - src/cosmology/cosmology_functions.cpp | 6 +----- src/global/global.h | 4 +--- 3 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index de2cf9d9d..e12ed39b1 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -21,7 +21,6 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part w0 = P->w0; wa = P->wa; - if (strcmp(P->init, "Read_Grid") == 0) { // Read scale factor value from Particles current_z = Particles.current_z; diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index dc8100ac4..6da05822c 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -5,9 +5,6 @@ #include "../grid/grid_enum.h" #include "../io/io.h" - - - void Grid3D::Initialize_Cosmology(struct Parameters *P) { chprintf("Initializing Cosmology... \n"); @@ -141,10 +138,9 @@ void Grid3D::Change_GAS_Frame_System(bool forward) } -// writing the expansion history - void Create_Expansion_History_File(struct Parameters P) { + // writing the expansion history if (not Is_Root_Proc()) { return; } diff --git a/src/global/global.h b/src/global/global.h index 3d73f43a3..d26d326eb 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -326,9 +326,7 @@ struct Parameters { Real wa; Real Init_redshift; Real End_redshift; - char scale_outputs_file[MAXLEN]; // File for the scale_factor output values - // for cosmological simulations - + char scale_outputs_file[MAXLEN]; // File for the scale_factor output values for cosmological simulations #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" #endif // COSMOLOGY #ifdef TILED_INITIAL_CONDITIONS From e4d824a06f119fdab293c63a5b18c725f19a38a1 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:14:18 -0700 Subject: [PATCH 08/62] Scope issues regarding expansion history log file functions. --- src/cosmology/cosmology_functions.cpp | 18 ++++++++++++++---- src/grid/grid3D.h | 2 +- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 6da05822c..69a83db1f 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -26,7 +26,7 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Real Cosmology::Get_da_from_dt(Real dt) { Real a2 = current_a * current_a; - Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real fac_de = pow(current_a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - current_a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return a_dot * dt; } @@ -34,7 +34,7 @@ Real Cosmology::Get_da_from_dt(Real dt) Real Cosmology::Get_dt_from_da(Real da) { Real a2 = current_a * current_a; - Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real fac_de = pow(current_a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - current_a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; return da / a_dot; } @@ -137,8 +137,18 @@ void Grid3D::Change_GAS_Frame_System(bool forward) } } +/* local function that designates whether we are using a root-process. It gives + * gives a sensible result regardless of whether we are using MPI */ +static inline bool Is_Root_Proc() +{ +#ifdef MPI_CHOLLA + return procID == root; +#else + return true; +#endif +} -void Create_Expansion_History_File(struct Parameters P) +void Cosmology::Create_Expansion_History_File(struct Parameters P) { // writing the expansion history if (not Is_Root_Proc()) { @@ -170,7 +180,7 @@ void Create_Expansion_History_File(struct Parameters P) out_file.close(); } -void Write_Expansion_History_Entry(void) +void Cosmology::Write_Expansion_History_Entry(void) { if (not Is_Root_Proc()) { return; diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index cb4c0dbbb..4a1fa17a5 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -854,7 +854,7 @@ class Grid3D void Change_DM_Frame_System(bool forward); void Change_GAS_Frame_System(bool forward); void Change_GAS_Frame_System_GPU(bool forward); - void Change_Cosmological_Frame_Sytem(bool forward); + void Change_Cosmological_Frame_System(bool forward); void Advance_Particles_KDK_Cosmo_Step1_function(part_int_t p_start, part_int_t p_end); void Advance_Particles_KDK_Cosmo_Step2_function(part_int_t p_start, part_int_t p_end); Real Calc_Particles_dt_Cosmo_function(part_int_t p_start, part_int_t p_end); From e2b5b16f40c21036bba28d5708b8e082b5767bb5 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:22:45 -0700 Subject: [PATCH 09/62] Sub to_string() for string() --- src/cosmology/cosmology.h | 4 ++-- src/cosmology/cosmology_functions.cpp | 6 +++--- src/gravity/gravity_functions.cpp | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cosmology/cosmology.h b/src/cosmology/cosmology.h index 53252ab16..020ea8c9b 100644 --- a/src/cosmology/cosmology.h +++ b/src/cosmology/cosmology.h @@ -71,8 +71,8 @@ class Cosmology Real Get_dt_from_da(Real da); // write expansion history log file - void Create_Expansion_History_File(struct Parameters P); - void Write_Expansion_History_Entry(Real t, Real a); + void Create_Expansion_History_File(struct Parameters *P); + void Write_Expansion_History_Entry(void); }; #endif diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 69a83db1f..f8efd03f9 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -170,8 +170,8 @@ void Cosmology::Create_Expansion_History_File(struct Parameters P) // convert now to string form char *dt = ctime(&now); - std::string message = std::string(Omega_M) + " " + std::string(Omega_K); - message += " " + std::string(Omega_R) + " " + std::string(Omega_L) + " " + std::string(w0) + " " + std::string(wa); + std::string message = std::to_string(Omega_M) + " " + std::to_string(Omega_K); + message += " " + std::to_string(Omega_R) + " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); @@ -186,7 +186,7 @@ void Cosmology::Write_Expansion_History_Entry(void) return; } - std::string message = std::string(Cosmo.t_secs / MYR) + " " + std::string(Cosmo.current_a); + std::string message = std::to_string(Cosmo.t_secs / MYR) + " " + std::to_string(Cosmo.current_a); std::string file_name(EXPANSION_HISTORY_FILE_NAME); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index aea3acb15..7e6cdb578 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -128,7 +128,7 @@ void Grid3D::set_dt_Gravity() Particles.dt = dt_physical; // Write expansion history - Write_Expansion_History_Entry(); + Cosmo.Write_Expansion_History_Entry(); #else // Not Cosmology // If NOT using COSMOLOGY From 66b59fa8609118d53ba3624d3de0ed89a71ceb52 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:25:41 -0700 Subject: [PATCH 10/62] Correct scope issues for cosmological expansion history log writting. --- src/cosmology/cosmology_functions.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index f8efd03f9..7d6931a88 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -11,7 +11,7 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) Cosmo.Initialize(P, Grav, Particles); // Create expansion history log file - Create_Log_File(P, Cosmo); + Cosmo.Create_Expansion_History_File(P); // Change to comoving Cosmological System Change_Cosmological_Frame_System(true); @@ -148,7 +148,7 @@ static inline bool Is_Root_Proc() #endif } -void Cosmology::Create_Expansion_History_File(struct Parameters P) +void Cosmology::Create_Expansion_History_File(struct Parameters *P) { // writing the expansion history if (not Is_Root_Proc()) { @@ -186,7 +186,7 @@ void Cosmology::Write_Expansion_History_Entry(void) return; } - std::string message = std::to_string(Cosmo.t_secs / MYR) + " " + std::to_string(Cosmo.current_a); + std::string message = std::to_string(t_secs / MYR) + " " + std::to_string(current_a); std::string file_name(EXPANSION_HISTORY_FILE_NAME); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); From 999e2f476657b16d28920e5cd2e2649807271b9e Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:28:24 -0700 Subject: [PATCH 11/62] Aggregate type error for std::ofstream out_file Not sure why this is triggered for cosmology_functions.cpp but not io.cpp. --- src/cosmology/cosmology_functions.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 7d6931a88..e2c2fd854 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -1,4 +1,7 @@ #ifdef COSMOLOGY + #include + #include + #include #include "../global/global.h" #include "../grid/grid3D.h" From 3a384bcffaab0550a1e9e96cbd3d30062c5584f4 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:30:27 -0700 Subject: [PATCH 12/62] Apparently I need to include fstream instead. --- src/cosmology/cosmology_functions.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index e2c2fd854..ac4052304 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -1,7 +1,5 @@ #ifdef COSMOLOGY - #include - #include - #include + #include #include "../global/global.h" #include "../grid/grid3D.h" From a0c3457a421f08a21375900be2bac643cbedc309 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:33:11 -0700 Subject: [PATCH 13/62] Typo in Advance_Particles_KDK_Step2_Cosmo_Kernel def extra comma deleted --- src/particles/particles_dynamics_gpu.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index fcbcdf335..594d2230f 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -237,7 +237,7 @@ __global__ void Advance_Particles_KDK_Step1_Cosmo_Kernel(part_int_t n_local, Rea __global__ void Advance_Particles_KDK_Step2_Cosmo_Kernel(part_int_t n_local, Real da, Real *vel_x_dev, Real *vel_y_dev, Real *vel_z_dev, Real *grav_x_dev, Real *grav_y_dev, Real *grav_z_dev, Real current_a, Real H0, Real cosmo_h, - Real Omega_M, Real Omega_L, Real Omega_K, Real, Omega_R, + Real Omega_M, Real Omega_L, Real Omega_K, Real Omega_R, Real w0, Real wa) { part_int_t tid = blockIdx.x * blockDim.x + threadIdx.x; From 858178471c64a43104e486801ffd15a433abff41 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:38:54 -0700 Subject: [PATCH 14/62] Adding missed line to expansion_history.txt file --- src/cosmology/cosmology_functions.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index ac4052304..4bfa9542b 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -178,6 +178,7 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) out_file.open(file_name.c_str(), std::ios::app); out_file << "\n"; out_file << "Run date: " << dt; + out_file << message.c_str() << std::endl; out_file.close(); } From f1966393a14d08ee479f90464089aec4ecb66834 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:44:16 -0700 Subject: [PATCH 15/62] Reordering expansion history header --- src/cosmology/cosmology_functions.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 4bfa9542b..0fc04c39f 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -171,8 +171,9 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) // convert now to string form char *dt = ctime(&now); - std::string message = std::to_string(Omega_M) + " " + std::to_string(Omega_K); - message += " " + std::to_string(Omega_R) + " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); + std::string message = std::to_string(H0) + " " + std::to_string(Omega_M); + message += " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); + message += " " + std::to_string(Omega_R) + " " + std::to_string(Omega_K); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); From 9fb92d64c6e87716c809c0ae6676ea22a1025e74 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:48:03 -0700 Subject: [PATCH 16/62] Print H0 in per mpc, not per kpc --- src/cosmology/cosmology_functions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 0fc04c39f..4ae4ed6dc 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -171,7 +171,7 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) // convert now to string form char *dt = ctime(&now); - std::string message = std::to_string(H0) + " " + std::to_string(Omega_M); + std::string message = std::to_string(H0 * 1e3) + " " + std::to_string(Omega_M); message += " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); message += " " + std::to_string(Omega_R) + " " + std::to_string(Omega_K); From e8147d8b73eaa0fbed2bb68f582f17af7189265a Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:51:08 -0700 Subject: [PATCH 17/62] Adding Omega_b to header for convenience --- src/cosmology/cosmology_functions.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 4ae4ed6dc..9f3724bbb 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -171,14 +171,15 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) // convert now to string form char *dt = ctime(&now); - std::string message = std::to_string(H0 * 1e3) + " " + std::to_string(Omega_M); + std::string message = "# H0 OmegaM Omega_b OmegaL w0 wa Omega_R Omega_K\n" + message += "# " + std::to_string(H0 * 1e3) + " " + std::to_string(Omega_M); + message += " " + std::to_string(Omega_b); message += " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); message += " " + std::to_string(Omega_R) + " " + std::to_string(Omega_K); std::ofstream out_file; out_file.open(file_name.c_str(), std::ios::app); - out_file << "\n"; - out_file << "Run date: " << dt; + out_file << "# Run date: " << dt; out_file << message.c_str() << std::endl; out_file.close(); } From 85ec02299c00ecca30a55a282bfc8e95be3c0a7b Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 16:53:36 -0700 Subject: [PATCH 18/62] Missing semicolon added --- src/cosmology/cosmology_functions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 9f3724bbb..fbbb2ce7e 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -171,7 +171,7 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) // convert now to string form char *dt = ctime(&now); - std::string message = "# H0 OmegaM Omega_b OmegaL w0 wa Omega_R Omega_K\n" + std::string message = "# H0 OmegaM Omega_b OmegaL w0 wa Omega_R Omega_K\n"; message += "# " + std::to_string(H0 * 1e3) + " " + std::to_string(Omega_M); message += " " + std::to_string(Omega_b); message += " " + std::to_string(Omega_L) + " " + std::to_string(w0) + " " + std::to_string(wa); From 80e0d0e25da1afe97a89c123e633c80691afcb47 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 22:41:28 -0700 Subject: [PATCH 19/62] Adding a more accurate cosmological time integrator. --- src/cosmology/cosmology.h | 2 ++ src/cosmology/cosmology_functions.cpp | 36 +++++++++++++++++++++++++++ src/gravity/gravity_functions.cpp | 5 ++-- 3 files changed, 41 insertions(+), 2 deletions(-) diff --git a/src/cosmology/cosmology.h b/src/cosmology/cosmology.h index 020ea8c9b..a82858788 100644 --- a/src/cosmology/cosmology.h +++ b/src/cosmology/cosmology.h @@ -67,6 +67,8 @@ class Cosmology Real Get_Hubble_Parameter(Real a); + Real dtda_cosmo(Real da, Real a); + Real Get_dt_from_da_rk(Real dt); Real Get_da_from_dt(Real dt); Real Get_dt_from_da(Real da); diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index fbbb2ce7e..c10624f80 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -24,6 +24,42 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) chprintf("Cosmology Successfully Initialized. \n\n"); } +Real Cosmology::dtda_cosmo(Real da, Real a) +{ + // Return dt/da * da + Real a2 = a * a; + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - current_a)); + Real a_dot = sqrt(Omega_R / a2 + Omega_M / a + a2 * Omega_L * fac_de + Omega_K) * H0; + return da / a_dot; +} + +Real Cosmology::Get_dt_from_da_rk(Real da) +{ + // Return dt/da * da + // But compute dt/da using a Runge-Kutta integration step + Real a3 = 0.3; + Real a4 = 0.6; + Real a5 = 1.0; + Real a6 = 0.875; + Real c1 = 37.0 / 378.0; + Real c3 = 250.0 / 621.0; + Real c4 = 125.0 / 594.0; + Real c6 = 512.0 / 1771.0; + + + // compute RK average derivatives + Real ak1 = dtda_cosmo(da, current_a); + Real ak3 = dtda_cosmo(da, current_a + a3 * da); + Real ak4 = dtda_cosmo(da, current_a + a4 * da); + Real ak6 = dtda_cosmo(da, current_a + a6 * da); + + // compute timestep + Real dt = (c1 * ak1 + c3 * ak3 + c4 * ak4 + c6 * ak6); + + // return timestep + return dt; +} + Real Cosmology::Get_da_from_dt(Real dt) { Real a2 = current_a * current_a; diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 7e6cdb578..4bf23dc72 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -54,7 +54,7 @@ void Grid3D::set_dt_Gravity() // Here da_min is the minumum between da_particles and da_hydro Real da_hydro; da_hydro = - Cosmo.Get_da_from_dt(dt_hydro) * Cosmo.current_a * Cosmo.current_a / Cosmo.H0; // Convet delta_t to delta_a + Cosmo.Get_da_from_dt(dt_hydro) * Cosmo.current_a * Cosmo.current_a / Cosmo.H0; // Convert delta_t to delta_a da_min = fmin(da_hydro, da_particles); // Find the minumum delta_a chprintf(" Delta_a_particles: %f Delta_a_gas: %f \n", da_particles, da_hydro); @@ -121,7 +121,8 @@ void Grid3D::set_dt_Gravity() #endif // Compute the physical time - dt_physical = Cosmo.Get_dt_from_da(Cosmo.delta_a); + // dt_physical = Cosmo.Get_dt_from_da(Cosmo.delta_a); + dt_physical = Cosmo.Get_dt_from_da_rk(Cosmo.delta_a); Cosmo.dt_secs = dt_physical * Cosmo.time_conversion; Cosmo.t_secs += Cosmo.dt_secs; chprintf(" t_physical: %f Myr dt_physical: %f Myr\n", Cosmo.t_secs / MYR, Cosmo.dt_secs / MYR); From 987181e5f5b10b1b92b7a7316287c80a9afcadc8 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 23:42:10 -0700 Subject: [PATCH 20/62] Updated time integrator for cosmological evolution. Added a computation of the offset in time at the start of the simulation, and a higher accuracy Runge-Kutta timestep scheme to improve the time integration given the typical scale factor step da. --- src/cosmology/cosmology.cpp | 17 ++++++++++++++++- src/cosmology/cosmology.h | 4 ++-- src/cosmology/cosmology_functions.cpp | 21 +++++++++++---------- src/gravity/gravity_functions.cpp | 7 +++---- src/particles/particles_dynamics.cpp | 4 ++-- 5 files changed, 34 insertions(+), 19 deletions(-) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index e12ed39b1..18577153f 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -45,8 +45,23 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part delta_a = max_delta_a; // Initialize Time and set the time conversion - t_secs = 0; time_conversion = KPC; + t_secs = 0; + + // get time offset + Real da_t_sec = 1.0e-2 * current_a; + Real a_t_sec = 1.0e-6; // small but non-zero + Real dt_physical; + while(a_t_sec < current_a){ + if(a_t_sec + da_t_sec > current_a){ + da_t_sec = current_a - a_t_sec; + } + dt_physical = Get_dt_from_da(da_t_sec, a_t_sec); + t_secs += dt_physical * time_conversion; + a_t_sec += da_t_sec; + // chprintf(" Revised a_t_sec %f da_t_sec %f t_start : %e Myr\n", a_t_sec, da_t_sec, t_secs / MYR); + } + chprintf(" Revised t_start : %f Myr\n", t_secs / MYR); // Set Normalization factors r_0_dm = P->xlen / P->nx; diff --git a/src/cosmology/cosmology.h b/src/cosmology/cosmology.h index a82858788..e941ebc23 100644 --- a/src/cosmology/cosmology.h +++ b/src/cosmology/cosmology.h @@ -68,9 +68,9 @@ class Cosmology Real Get_Hubble_Parameter(Real a); Real dtda_cosmo(Real da, Real a); - Real Get_dt_from_da_rk(Real dt); + Real Get_dt_from_da_rk(Real da, Real a); Real Get_da_from_dt(Real dt); - Real Get_dt_from_da(Real da); + Real Get_dt_from_da(Real da, Real a); // write expansion history log file void Create_Expansion_History_File(struct Parameters *P); diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index c10624f80..0ca3eaa91 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -33,7 +33,7 @@ Real Cosmology::dtda_cosmo(Real da, Real a) return da / a_dot; } -Real Cosmology::Get_dt_from_da_rk(Real da) +Real Cosmology::Get_dt_from_da_rk(Real da, Real a) { // Return dt/da * da // But compute dt/da using a Runge-Kutta integration step @@ -48,10 +48,10 @@ Real Cosmology::Get_dt_from_da_rk(Real da) // compute RK average derivatives - Real ak1 = dtda_cosmo(da, current_a); - Real ak3 = dtda_cosmo(da, current_a + a3 * da); - Real ak4 = dtda_cosmo(da, current_a + a4 * da); - Real ak6 = dtda_cosmo(da, current_a + a6 * da); + Real ak1 = dtda_cosmo(da, a); + Real ak3 = dtda_cosmo(da, a + a3 * da); + Real ak4 = dtda_cosmo(da, a + a4 * da); + Real ak6 = dtda_cosmo(da, a + a6 * da); // compute timestep Real dt = (c1 * ak1 + c3 * ak3 + c4 * ak4 + c6 * ak6); @@ -68,12 +68,13 @@ Real Cosmology::Get_da_from_dt(Real dt) return a_dot * dt; } -Real Cosmology::Get_dt_from_da(Real da) +Real Cosmology::Get_dt_from_da(Real da, Real a) { - Real a2 = current_a * current_a; - Real fac_de = pow(current_a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - current_a)); - Real a_dot = sqrt(Omega_R / a2 + Omega_M / current_a + a2 * Omega_L * fac_de + Omega_K) * H0; - return da / a_dot; + return Get_dt_from_da_rk(da,a); + /* Real a2 = a * a; + Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); + Real a_dot = sqrt(Omega_R / a2 + Omega_M / a + a2 * Omega_L * fac_de + Omega_K) * H0; + return da / a_dot; */ } Real Cosmology::Get_Hubble_Parameter(Real a) diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 4bf23dc72..6edbc47f4 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -107,7 +107,7 @@ void Grid3D::set_dt_Gravity() // Set delta_a after it has been computed Cosmo.delta_a = da_min; // Convert delta_a back to delta_t - dt_min = Cosmo.Get_dt_from_da(Cosmo.delta_a) * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a); + dt_min = Cosmo.Get_dt_from_da(Cosmo.delta_a, Cosmo.current_a) * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a); // Set the new delta_t for the hydro step H.dt = dt_min; chprintf(" Current_a: %f delta_a: %f dt: %f\n", Cosmo.current_a, Cosmo.delta_a, H.dt); @@ -115,14 +115,13 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif // Compute the physical time - // dt_physical = Cosmo.Get_dt_from_da(Cosmo.delta_a); - dt_physical = Cosmo.Get_dt_from_da_rk(Cosmo.delta_a); + dt_physical = Cosmo.Get_dt_from_da(Cosmo.delta_a, Cosmo.current_a); Cosmo.dt_secs = dt_physical * Cosmo.time_conversion; Cosmo.t_secs += Cosmo.dt_secs; chprintf(" t_physical: %f Myr dt_physical: %f Myr\n", Cosmo.t_secs / MYR, Cosmo.dt_secs / MYR); diff --git a/src/particles/particles_dynamics.cpp b/src/particles/particles_dynamics.cpp index 43e482023..4a609a09e 100644 --- a/src/particles/particles_dynamics.cpp +++ b/src/particles/particles_dynamics.cpp @@ -81,7 +81,7 @@ Real Grid3D::Calc_Particles_dt_GPU() scale_factor = 1 / (Cosmo.current_a * Cosmo.Get_Hubble_Parameter(Cosmo.current_a)) * Cosmo.cosmo_h; vel_factor = Cosmo.current_a / scale_factor; da_min = vel_factor / max_dti; - dt_min = Cosmo.Get_dt_from_da(da_min); + dt_min = Cosmo.Get_dt_from_da(da_min, Cosmo.current_a); #else dt_min = 1 / max_dti; #endif @@ -377,7 +377,7 @@ Real Grid3D::Calc_Particles_dt_Cosmo_function(part_int_t p_start, part_int_t p_e da_min = fmin(Particles.G.dx / vx_max, Particles.G.dy / vy_max); da_min = fmin(Particles.G.dz / vz_max, da_min); da_min *= vel_factor; - dt_min = Cosmo.Get_dt_from_da(da_min); + dt_min = Cosmo.Get_dt_from_da(da_min, Cosmo.current_a); return Particles.C_cfl * dt_min; } From addf760b501ba5a879167adeb3823cb371159327 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 23:48:27 -0700 Subject: [PATCH 21/62] Clang formatting take 0 --- src/cosmology/cosmology.cpp | 9 +++++---- src/cosmology/cosmology_functions.cpp | 20 ++++++++++---------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index 18577153f..6643e8805 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -50,14 +50,15 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part // get time offset Real da_t_sec = 1.0e-2 * current_a; - Real a_t_sec = 1.0e-6; // small but non-zero + // small but non-zero + Real a_t_sec = 1.0e-6; Real dt_physical; - while(a_t_sec < current_a){ - if(a_t_sec + da_t_sec > current_a){ + while (a_t_sec < current_a) { + if (a_t_sec + da_t_sec > current_a) { da_t_sec = current_a - a_t_sec; } dt_physical = Get_dt_from_da(da_t_sec, a_t_sec); - t_secs += dt_physical * time_conversion; + t_secs += dt_physical * time_conversion; a_t_sec += da_t_sec; // chprintf(" Revised a_t_sec %f da_t_sec %f t_start : %e Myr\n", a_t_sec, da_t_sec, t_secs / MYR); } diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 0ca3eaa91..0016dd546 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -37,14 +37,14 @@ Real Cosmology::Get_dt_from_da_rk(Real da, Real a) { // Return dt/da * da // But compute dt/da using a Runge-Kutta integration step - Real a3 = 0.3; - Real a4 = 0.6; - Real a5 = 1.0; - Real a6 = 0.875; - Real c1 = 37.0 / 378.0; - Real c3 = 250.0 / 621.0; - Real c4 = 125.0 / 594.0; - Real c6 = 512.0 / 1771.0; + Real a3 = 0.3; + Real a4 = 0.6; + Real a5 = 1.0; + Real a6 = 0.875; + Real c1 = 37.0 / 378.0; + Real c3 = 250.0 / 621.0; + Real c4 = 125.0 / 594.0; + Real c6 = 512.0 / 1771.0; // compute RK average derivatives @@ -54,7 +54,7 @@ Real Cosmology::Get_dt_from_da_rk(Real da, Real a) Real ak6 = dtda_cosmo(da, a + a6 * da); // compute timestep - Real dt = (c1 * ak1 + c3 * ak3 + c4 * ak4 + c6 * ak6); + Real dt = (c1 * ak1 + c3 * ak3 + c4 * ak4 + c6 * ak6); // return timestep return dt; @@ -70,7 +70,7 @@ Real Cosmology::Get_da_from_dt(Real dt) Real Cosmology::Get_dt_from_da(Real da, Real a) { - return Get_dt_from_da_rk(da,a); + return Get_dt_from_da_rk(da, a); /* Real a2 = a * a; Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / a + a2 * Omega_L * fac_de + Omega_K) * H0; From 05893a2ca7bb5607b01c6ad0c4f62ea79ce95f8f Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 20 Jun 2024 23:58:27 -0700 Subject: [PATCH 22/62] Clang tidy take 1 --- src/cosmology/cosmology.cpp | 4 ++-- src/cosmology/cosmology_functions.cpp | 1 - src/global/global.cpp | 9 ++++----- src/global/global.h | 7 +++++-- src/gravity/gravity_functions.cpp | 4 ++-- 5 files changed, 13 insertions(+), 12 deletions(-) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index 6643e8805..06a69f794 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -51,13 +51,13 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part // get time offset Real da_t_sec = 1.0e-2 * current_a; // small but non-zero - Real a_t_sec = 1.0e-6; + Real a_t_sec = 1.0e-6; Real dt_physical; while (a_t_sec < current_a) { if (a_t_sec + da_t_sec > current_a) { da_t_sec = current_a - a_t_sec; } - dt_physical = Get_dt_from_da(da_t_sec, a_t_sec); + dt_physical = Get_dt_from_da(da_t_sec, a_t_sec); t_secs += dt_physical * time_conversion; a_t_sec += da_t_sec; // chprintf(" Revised a_t_sec %f da_t_sec %f t_start : %e Myr\n", a_t_sec, da_t_sec, t_secs / MYR); diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 0016dd546..e72ff2c5a 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -46,7 +46,6 @@ Real Cosmology::Get_dt_from_da_rk(Real da, Real a) Real c4 = 125.0 / 594.0; Real c6 = 512.0 / 1771.0; - // compute RK average derivatives Real ak1 = dtda_cosmo(da, a); Real ak3 = dtda_cosmo(da, a + a3 * da); diff --git a/src/global/global.cpp b/src/global/global.cpp index 01acc5542..3b95cf733 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,11 +102,10 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) -const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", - "w0", "wa", "Init_redshift", - "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; +const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", + "w0", "wa", "Init_redshift", + "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/global/global.h b/src/global/global.h index d26d326eb..c5d458562 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -326,9 +326,12 @@ struct Parameters { Real wa; Real Init_redshift; Real End_redshift; - char scale_outputs_file[MAXLEN]; // File for the scale_factor output values for cosmological simulations + + // File for the scale_factor output values for cosmological simulations + char scale_outputs_file[MAXLEN]; + #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" -#endif // COSMOLOGY +#endif #ifdef TILED_INITIAL_CONDITIONS Real tile_length; #endif // TILED_INITIAL_CONDITIONS diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 6edbc47f4..76dec80da 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,8 +115,8 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / - SLOW_FACTOR; + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / + (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif From 5d79a80d91f17fe677c84769ddb29aac87443c68 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:03:43 -0700 Subject: [PATCH 23/62] Clang formatting take 2 --- src/global/global.cpp | 9 +++++---- src/global/global.h | 2 +- src/gravity/gravity_functions.cpp | 5 +++-- src/system_tests/system_tester.cpp | 3 ++- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 3b95cf733..4594fad4f 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,10 +102,11 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) -const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", - "w0", "wa", "Init_redshift", - "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; +const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", + "phi", "theta", "delta", "nzr", "nxr", "H0", + "Omega_M", "Omega_L", "Omega_R", "w0", "wa", + "Init_redshift", "End_redshift", "tile_length", + "n_proc_x", "n_proc_y", "n_proc_z" }; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/global/global.h b/src/global/global.h index c5d458562..2a352e280 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -328,7 +328,7 @@ struct Parameters { Real End_redshift; // File for the scale_factor output values for cosmological simulations - char scale_outputs_file[MAXLEN]; + char scale_outputs_file[MAXLEN]; #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" #endif diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 76dec80da..a9ac7d9da 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,8 +115,9 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / - (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * + Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; + H.min_dt_slow = min_dt_slow; #endif diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 0abebdeb8..93292b5e8 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -355,7 +355,8 @@ void system_test::SystemTestRunner::launchCholla() } #ifdef COSMOLOGY try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", + _outputDirectory + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } From 0da6407a0021c9438669e7b3faac932344c42365 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:09:20 -0700 Subject: [PATCH 24/62] Clang formatting take 3 --- src/global/global.h | 5 ++--- src/gravity/gravity_functions.cpp | 4 +--- src/system_tests/system_tester.cpp | 4 ++-- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/global/global.h b/src/global/global.h index 2a352e280..fa7ee1b29 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -329,9 +329,8 @@ struct Parameters { // File for the scale_factor output values for cosmological simulations char scale_outputs_file[MAXLEN]; - -#define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" -#endif + #define EXPANSION_HISTORY_FILE_NAME "expansion_history.txt" +#endif // COSMOLOGY #ifdef TILED_INITIAL_CONDITIONS Real tile_length; #endif // TILED_INITIAL_CONDITIONS diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index a9ac7d9da..30054ef4e 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,9 +115,7 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * - Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; - + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 93292b5e8..d6a888c43 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -353,14 +353,14 @@ void system_test::SystemTestRunner::launchCholla() } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } -#ifdef COSMOLOGY + #ifdef COSMOLOGY try { std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } -#endif + #endif } // ============================================================================= From 45630db58a0f1b0b8af200b3e007f06aa088f4d4 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:15:38 -0700 Subject: [PATCH 25/62] Clang formatting take 4 --- src/global/global.cpp | 9 ++++----- src/gravity/gravity_functions.cpp | 3 ++- src/particles/particles_dynamics_gpu.cu | 4 ++-- src/system_tests/system_tester.cpp | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 4594fad4f..3af5a4ced 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,11 +102,10 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) -const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", - "phi", "theta", "delta", "nzr", "nxr", "H0", - "Omega_M", "Omega_L", "Omega_R", "w0", "wa", - "Init_redshift", "End_redshift", "tile_length", - "n_proc_x", "n_proc_y", "n_proc_z" }; +const std::set optionalParams = {"flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "w0", + "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", + "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 30054ef4e..f4194dbef 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,7 +115,8 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * + Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index 594d2230f..75a6028ed 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -18,8 +18,8 @@ // FUTURE FIX: The Hubble function was defined here because I couldn't get it // form other file, tried -dc flag when compiling buu paris broke. -__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, - Real Omega_K, Real Omega_R, Real w0, Real wa) +__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, + Real Omega_R, Real w0, Real wa) { Real a2 = a * a; Real a3 = a2 * a; diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index d6a888c43..826dac3c3 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -355,8 +355,8 @@ void system_test::SystemTestRunner::launchCholla() } #ifdef COSMOLOGY try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", - _outputDirectory + "/expansion_history.txt"); + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } From 4a4b15438a5da3f0849c14a51d5385bf05b19204 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:22:47 -0700 Subject: [PATCH 26/62] Clang formatting take 5 --- src/gravity/gravity_functions.cpp | 4 ++-- src/system_tests/system_tester.cpp | 7 +++---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index f4194dbef..f66cfd553 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,8 +115,8 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * - Cosmo.current_a) / SLOW_FACTOR; + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / + SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 826dac3c3..474a1f398 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -353,14 +353,13 @@ void system_test::SystemTestRunner::launchCholla() } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } - #ifdef COSMOLOGY +#ifdef COSMOLOGY try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + - "/expansion_history.txt"); + std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } - #endif +#endif } // ============================================================================= From dc2b517b995be954dfe636473007155579135471 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:26:09 -0700 Subject: [PATCH 27/62] Clang formatting take 6 Removing functionality just to get this to pass --- src/system_tests/system_tester.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/system_tests/system_tester.cpp b/src/system_tests/system_tester.cpp index 474a1f398..1888fd752 100644 --- a/src/system_tests/system_tester.cpp +++ b/src/system_tests/system_tester.cpp @@ -353,13 +353,6 @@ void system_test::SystemTestRunner::launchCholla() } catch (const std::filesystem::filesystem_error &error) { // This file might not exist and isn't required so don't worry if it doesn't exist } -#ifdef COSMOLOGY - try { - std::filesystem::rename(::globalChollaRoot.getString() + "/expansion_history.txt", _outputDirectory + "/expansion_history.txt"); - } catch (const std::filesystem::filesystem_error &error) { - // This file might not exist and isn't required so don't worry if it doesn't exist - } -#endif } // ============================================================================= From 499b5009fc40246aab944ae2a99778d91a9ec094 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:29:10 -0700 Subject: [PATCH 28/62] Clang formatting take 7 --- src/gravity/gravity_functions.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index f66cfd553..41fb6c44e 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,8 +115,9 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a) / Particles.C_cfl * Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / - SLOW_FACTOR; + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a); + min_dt_slow /= Particles.C_cfl; + min_dt_slow *= Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; #endif From 84fbe8b2ce85d50f6324f8b8650539d0ee4e771f Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:35:32 -0700 Subject: [PATCH 29/62] Clang formatting take 8 --- src/global/global.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 3af5a4ced..a209676ee 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,10 +102,11 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) -const std::set optionalParams = {"flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "w0", - "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", - "n_proc_y", "n_proc_z"}; +const std::set optionalParams = { + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", + "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", + "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From bedf95dd6c00b1eaaf7dc9502f632eaa95878794 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:39:09 -0700 Subject: [PATCH 30/62] Clang formatting 9 --- src/global/global.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index a209676ee..ffc2b3a6b 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,10 +103,10 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", - "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", - "n_proc_y", "n_proc_z"}; + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", + "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", + "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From e6a38754e2bcf052f322273791a9b60257e4237e Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:42:49 -0700 Subject: [PATCH 31/62] Clang formatting take 10 --- src/global/global.cpp | 2 +- src/gravity/gravity_functions.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index ffc2b3a6b..8050810fc 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,7 +103,7 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", + "flag_delta", "ddelta_dt","n_delta","Lz","Lx","phi","theta", "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; diff --git a/src/gravity/gravity_functions.cpp b/src/gravity/gravity_functions.cpp index 41fb6c44e..7f1f9cc28 100644 --- a/src/gravity/gravity_functions.cpp +++ b/src/gravity/gravity_functions.cpp @@ -115,7 +115,7 @@ void Grid3D::set_dt_Gravity() #ifdef AVERAGE_SLOW_CELLS // Set the min_delta_t for averaging a slow cell da_particles = fmin(da_particles, Cosmo.max_delta_a); - min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a); + min_dt_slow = Cosmo.Get_dt_from_da(da_particles, Cosmo.current_a); min_dt_slow /= Particles.C_cfl; min_dt_slow *= Cosmo.H0 / (Cosmo.current_a * Cosmo.current_a) / SLOW_FACTOR; H.min_dt_slow = min_dt_slow; From 5c959b95bce72f5a78031f85522d011f35eb55a2 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:45:24 -0700 Subject: [PATCH 32/62] Clang formatting take 11 --- src/global/global.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 8050810fc..522e34eea 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,10 +103,9 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { - "flag_delta", "ddelta_dt","n_delta","Lz","Lx","phi","theta", - "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", - "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", - "n_proc_y", "n_proc_z"}; + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", + "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", + "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From e030afb6a40cb0b3a72fca6ffcc4a0d026695c2e Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:50:00 -0700 Subject: [PATCH 33/62] Clang formatting take 12 --- src/global/global.cpp | 6 +++--- src/particles/particles_dynamics_gpu.cu | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 522e34eea..61f87e23c 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,9 +103,9 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", - "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", - "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", + "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", + "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/particles/particles_dynamics_gpu.cu b/src/particles/particles_dynamics_gpu.cu index 75a6028ed..854541ac3 100644 --- a/src/particles/particles_dynamics_gpu.cu +++ b/src/particles/particles_dynamics_gpu.cu @@ -18,8 +18,8 @@ // FUTURE FIX: The Hubble function was defined here because I couldn't get it // form other file, tried -dc flag when compiling buu paris broke. -__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, - Real Omega_R, Real w0, Real wa) +__device__ Real Get_Hubble_Parameter_dev(Real a, Real H0, Real Omega_M, Real Omega_L, Real Omega_K, Real Omega_R, + Real w0, Real wa) { Real a2 = a * a; Real a3 = a2 * a; From fbf9cca1ef0e5ea7a83b08490d5571733e6e06dd Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:52:26 -0700 Subject: [PATCH 34/62] Clang formatting take 13 --- src/cosmology/cosmology_functions.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index e72ff2c5a..7a1ee8afe 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -178,11 +178,11 @@ void Grid3D::Change_GAS_Frame_System(bool forward) * gives a sensible result regardless of whether we are using MPI */ static inline bool Is_Root_Proc() { -#ifdef MPI_CHOLLA + #ifdef MPI_CHOLLA return procID == root; -#else + #else return true; -#endif + #endif } void Cosmology::Create_Expansion_History_File(struct Parameters *P) From e3ac9f744bad34416d04dc95810f42c6f182d6b9 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 00:59:08 -0700 Subject: [PATCH 35/62] Clang formatting take 14 NOLINTBEGIN AND NOLINTEND --- src/global/global.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/global/global.cpp b/src/global/global.cpp index 61f87e23c..294e7f5dc 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,10 +102,12 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) +// NOLINTBEGIN(*) const std::set optionalParams = { "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; +// NOLINTEND(*) bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From 5eaf89eb9e8edc54b8144424aac264daca5ae6bb Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 01:00:54 -0700 Subject: [PATCH 36/62] Clang formatting take 16 --- src/global/global.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 294e7f5dc..76abcc3a8 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,10 +103,10 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) // NOLINTBEGIN(*) -const std::set optionalParams = { - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", - "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", - "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; +const std::set optionalParams = { // NOLINT + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", // NOLINT + "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", // NOLINT + "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; // NOLINT // NOLINTEND(*) bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From ee7e89044aee9e76c7c66fdef8186a93922772ba Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 01:04:10 -0700 Subject: [PATCH 37/62] Clang formatting take 17 --- src/global/global.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/global/global.cpp b/src/global/global.cpp index 76abcc3a8..aff1bfa55 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -102,12 +102,7 @@ char *Trim(char *s) } // NOLINTNEXTLINE(cert-err58-cpp) -// NOLINTBEGIN(*) -const std::set optionalParams = { // NOLINT - "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", "nzr", "nxr", // NOLINT - "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", "wa", "Init_redshift", // NOLINT - "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; // NOLINT -// NOLINTEND(*) +const std::set optionalParams = { "flag_delta","ddelta_dt","n_delta","Lz","Lx","phi","theta","delta","nzr","nxr","H0","Omega_M","Omega_L","Omega_R","Omega_K","w0","wa","Init_redshift","End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z" }; // NOLINT bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From 81207de8a3016875469aa431509348003086c9b7 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 01:06:42 -0700 Subject: [PATCH 38/62] Clang formatting take 18 --- src/global/.global.cpp.swp | Bin 0 -> 28672 bytes src/global/global.cpp | 1 + 2 files changed, 1 insertion(+) create mode 100644 src/global/.global.cpp.swp diff --git a/src/global/.global.cpp.swp b/src/global/.global.cpp.swp new file mode 100644 index 0000000000000000000000000000000000000000..80b64b8ba50f1d0cb3d76c0c5a98ae06f00f030b GIT binary patch literal 28672 zcmeI432-D=d4LBD1QW0^hT=4SmhJUOmNZ(M)81WME6uKjby?EdUUr?PHPbWG-k$DB zcaNmeVkjV$RJa`CB8Cu?0z0Hq6i1~vQ=!7)a0CK5On^WrAQYr9l}d3Ygh0Okb1s1AySZl5qN;S`M%tGLo3RTBjH=X#;WU0~EUC!|C zPy<5^oQ(#akh}V-{kp4PI&z7+@JYwcM#wPEPy<5^3^g#+z)%B24Gc9f)WA>!Lk;{N z)j-&IbnfNk{gHA0Ulsp-#E!p@ihtiUAl#$l>qh+df%tD*&ST==rTFjF@s0NVbK>jQ z$A5os$Nl*7@Xt^KLk$cyFx0?M149iAH89k`Py<5^3^g#+z)%B24V;Ar%yKSw5%H$Q z0U-PTv5W`L&*eS_ABDHW+u-%^Qn29|Ov9Bh2G4|Nz}L>peVt6+E^zpgeU&1@!K6oR%5)Q-la2;F%Km65P?w{bR@Co=MXhHx! zti$E-6u1x`1wVcq>BHCH%kU5IMc9BF;YBb3KRcJ}@BsKQ4F};Va3MSb{`9fA+&keO zxD#FhFNc>w3G(n{coIDLE4kdq;bU+Qgm61lUVv7 z4uftNJKNQzmRq(fIi(t=pG*Fz2NQ<^&Vg8l*9;rYFfb~%W1X0~O&wHMoSeGiTH+K6 zYRHjx$AX4-#vs@szfu+KlCHan`1$ zEvr(!JVAW95qesrqCYf_CBa)q%= z)X9pg&N7dp%6dR4qe4paipprL>ST-stG;geL0%WU-CXwRl8?t(=`Xa%TlYLiEjeWb zwP@AEZ>uj~^4uV#=0%6nfX7wx%7*DQximn6j$E+G<%L&bq@sYCrVw6L^yPd1^+x(p0;=e}C=6~9H6grOOZu=5`Oa8xM zivO1CR@F2XI^X6y-&Q)`Zsgmh1m2W{wQ_x{9MtTJBwlpOM)#*032EboW4YC^MuR4A z4c{vnWU2GLmH56rqVBj>hJ?i>+WT6r^KG$s>zed{{Mcy9#ggTR6PE8k_XQK278}jw zq<==eQw7pL(R-=U%F}RzsE-!L)Qs0?(FR>RsL=}vA^nubNARe72m7zQQX+5GWKhcI z?I38FB`ZIMj8V03pUMaO)d3Y;#V_f0$XxU7x^k>ar~cuHs zH8GK%g2>bo)452V6fGVif}W0k)CvPKQC40>g+61N9`BIi zFapST^cPxMbv>P1i|v&f)tl#qN|EO&$)hXnMSksR8auQa4#eDJ0iSc+3}imY}J&r%49P<~pN<);SYcbYfB=e8RU7 zO@s~8x6O6O3Rs{qSo>)8TGUz?;}qv9i6>XILaiVm;}S?Kn^1)@w<#l}1QV4L!qQV? zoDQ0+$5Emhq0f%R0hP$-rgT%8TN6cUXK6>-bdTPx1S*f97W0jTI!X?I%c^GaYuzfD z&1epfZg*8wR+h{$Qqj?Hkesrk#&n0W93r2lU3Er_3Ej+OQr)6!Pf4#ktyQyHb@NPV zj@s&J{0k;ot5cV(a;VoTvs&ErCYb?cc@RXG(d1T7^A~+1v}dde?YiYph8sp;m1M-2 zoU$gKCoQUdox+x+A5^VSu6dQR*~;(Jx5uLefZUB|8M(~76UCM(zw9+>5Q+|KtZ^?d z(`!#@+BuebD>Z%JB9E-EN~g4pp6sW~(j7SNkJGEG#g(;u)vTN9kd%cX|EO8#SG~zh zNq4NW?qQgb4&y4%wI-zlQaFimQ2x|S!`AC7+LE+(eTjnTb-^md(DH()V%p(*N<^S2 z7X5)4n32Kr)rw7(xeN}TvqADS)}95&V8PwRZ9Mo$B`zAX(`(bla&g62o?cm-n_#!GFt42*DQz*2fwHM=k5%pp!Uj21NP{=4| z*D5eDBy!#82WX0idVaTpf^@dY+M#49y`tE8jgqbRU6UdtqcWh>dasniRzu70CLn5y zM$Ps~`PvAJm+lHY)_QwDii*V)BGZm;^Q_giofApX2kpiFKL=a#nb?|Q|6d>5dY{9V zzZYHuCtw^dg=fKekb}>km&^SzycO<-*TObzKou+~!w6gq-@>;4I=lhyghTLj_z|}J z7vLnEfE(b6a4tL+K8fAF4J{ah?_-}o2ycbk;c4(x_ysokU&9Hw1ulTc!H==OKMgO2 zo8f720sI$b{#z)+OW-oN2t*gpfZ=DTfuRP58W?Ke|D^`*=&t@)PPJW|rf*ank44%+ z?Hruc%S^-b^CRN;Y~RJjF*P}*SGSQJ(q5G0M{b#3S)5xudQEq^5jn=no)x%PgskkM zRZ>!1H|(UEwJK)Q(N0D@i8^^%giT$72T>;|#*A2!;trjxnSnv>&H9e~Y5zjJ{9B!w zo}XSZj?6DDt@M?-Y`K9Qwgwft=k~);>K=IyOR;+-jwH)HTU=b7TRT3W;32MP%QwYG zIjHb`K^_JL^vQf!l+Y)&M-gkq1+1CrwPP#ASQ%Igy;fb^1GLT6vF@T@x_Qa6`b>6S zkBozu%chgpjZ`<{1CfP!C_kdbl0mTzS)I2%^i8RXaaA{4n39`T+tq^`VncjyrmT0R zXLS3-a|>5%drAuwKVO^!*9bzZ5yYDm{XMz&zErBtOB|c64T;`1Ka%4>w+wM0dZaFc zx$gN`Z*YE$=3QF9JwbLF-)}}FtjiW2mH}{uNwk2d3b+>c@N3}WJx!6j(xE)nHl^Uqk8)z#v&{hlXc zjoei|EPK)tH(A&;9j7H;rH&()5S~+>2xOn5*W=|kE$P=i-noj2Urt|bd+0i?vGb+Vd1;ZhYxuaDo3b%t zhMt|*14uO6MiOKWls$lYBOmR}j3lFt;F+g~MK1QZ<<1^kUBlMa8;;S0*qsU!nLwy3 zI}34B4o+E{_*}}iKOPM`<7#30*7@S%?wRREOUugAv9;x6YsTu0#aoIit0=2ecBf?; z%cqR8;doW!R8;!PUP|8^F_Ea){}-^ltIK(N?Em)p{tsc_zX>h_vHO37{r*FE1Bh*3 zf#<{X;B(mUe+8d`JK!34KKvWD`}g5@;del6`xn6Z@NMk)`{8boJ^Lak_*d-sPr>Wp z*I+-K4`0KMe;=HN1rXc*8`$voK>*K($H3>Y+usIv!$}y2XTyW|I6eU%f|tULa6WvU zee<`&ufsYN;K}e7_RsV16ZXhI0v;^EQ{m_Ai~kAy2HXUr@Jsaf4fp`u4Wds6dOkG! zLCtf!ZnJz#G|Nt%x44s`rX_f~psj3zBzU@@_VgHSN9iuOr$*QlFI}|5q*ju}rWwlY zOuwS&VZ9-P!_Ema&-GYgjH>Nl`$w59QpB`<(EAmErzWG(Td}wZMSg~Xq+fP=QOG_g zAPm#lG+Q!IWiLgqgaHL$0*dxdJ3FR0*~{iS+j3!Rf}N*^V^>-evXRZ+)P%Lw@Nt;) zA>v%s^?Yl=Ecy6x^^S9UQN;?2t%ow|VX0%5C)+W~O%70aIMO+=v1!?UndO+oMl3FD zy#W{EhMr``VNa}7v)D?&1i|^w-b~pJWP>W+%MnMn>uoZHh|R;cgXvZ$Lc82j`(Gf_ zX*Ov-1}MlxlcVV!XHcF|K?i6#+E z9$McEH_U3`@#H}au`JB+>Z|ZLS6vH#zai_Ez0PUr0upH}KD#9*>z{{|7)zzJF)KA= zvW|L4iLvDBBo@!q{Tl7*a~YEB#!IH64?5|(n=KBix}#=LbS=Mnw#w|pPbM5KL`#U6;WZ2ocV^kiI_ zCEJ5Yh6H-z4H6wVcFD>lRL2-3&i6t&s2~S|ZnvV-8halQSn>BoAC%BpCeV#GP;k5@ zn@J~(zj_~ckkDGI%mkW66FSBKae99*gD?rdL(~ovM|PUijGKNFOG5U@VkcN1y3!MK z$zy+;!-b3}eoJP26w^M086RL^lk)Wc+>yJEHGDEUDJFa4EN3W(tP?u9QJZoSm9vT3 zM3z<1BqulJ;D|>!f71>umP2Fw>l}St?z4-<#e*HKg{pN{E}WL)JP_s8*0Z)(m9|`?!&o}PE|xyW3B{roZ&@fFGC6HQ*5w=$YohpMocP*A@lxE1_NK9r;}Z<@YqG+j|r|h0qnkf>8VwWTxw~7}r^p+pe`KaYWG_cLir3#Ot7|BRT-e;oDQu zk>3h$oAWkf&zsJC!wdpLTBB_RM(A;9QQICO)QWxP6*tPf*(4X7Pd%tk?7M9*CqFrJ zcSz~Oqo>>2Pd`N(c_-rPC22Xv%8|&ZE(hvN*5f3_ro%|1!K6?t)XW2>ak%_!f5mCqef4UjfJA06ZCfiv9mLAU=WL z1{%oGfa1s0hU%&(Kes~QyFb$W$Pw@wQ8Qu+YzTSc(a1s0xpTHO4 z1Mpg?!wg&o|AC+2EAS3@85H5E@O69z?}bg6gQvj{@D+R%-V9Z^9?paBqW?b!sfV{n zUGPzSbd#7k)UZSRYO zFNZSLGV#(ymu<6=rSs_U>7ovA8{(0|uG~}RB=`V9t!#n_K0we`HbDeW7qpy74x!V< zY-H1r1Ro$Mn}#+9YG@;yh9r1;j%Z++bfm+li=tsT+1imaKu{~2Ac7ANw3ST|!P5mT zXOctcbTKrHEEsNZFU&@kbcjFTAb|et5$G4k?2#xf* zr~T(R#sG0nW^qxBbaAcBxr<{A5SKZ3t=)3Bl{t5DjP%@XWzJn3BVF7!&&*}$E)qOl z5a%ee2$JCGg0`{=lHlorWRH1g(Zv48oX)*lY){xXH#Ff0jKB}D`Q_aId*H25fdc#-d;bA=JDi49cpjV!-^JelGk6o6f;o5= z{0v+F@8JFLI%vQNxE3A<-^SkmDEtMy6K;nC@N{?(oBz}B`>+WOcp*Fi{vF%@es}}C z8UnZu9uJR(AL9dfAA~Ri7s6Mm6RDSX#r}XEkAR%$h>wwQnvKUXvv}(JBkyQ2JBDYI zN0D2Bo%x(gBHrE-F#6vfkZIsxLPK#c8klS%Yuz_^Bi;Ul0k|C4rx1FRmr1>nSJ8w%# z7m~&5=4~sfLU;j71}(|!T2jR@Rq%lL-p5HDK3$Z4V@w87(W_=sbH#gPveX!Fok?bnVj8#1EzTN8 zipAN(EbB;nmrS@kO2+Yv?72up9F&dK#cn2IRdq6zn3%gy{wKP5!*{E^k%BV|y4~$* z=;*{rWGGHIHH}Skqy3to3{5f-E0Inj>&wQvWAd_zy?2;6M!LA optionalParams = { "flag_delta","ddelta_dt","n_delta","Lz","Lx","phi","theta","delta","nzr","nxr","H0","Omega_M","Omega_L","Omega_R","Omega_K","w0","wa","Init_redshift","End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z" }; // NOLINT bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); From 7f885a945c27bb2ddc3d86f7fb70762cf727b7d6 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 21 Jun 2024 01:10:28 -0700 Subject: [PATCH 39/62] prevent .swp from vi --- .gitignore | 1 + src/global/.global.cpp.swp | Bin 28672 -> 0 bytes 2 files changed, 1 insertion(+) delete mode 100644 src/global/.global.cpp.swp diff --git a/.gitignore b/.gitignore index 864a8ab2c..87e13814b 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ bin/* *.dll *.exe *.o +*.swp *.so *.mod *.a diff --git a/src/global/.global.cpp.swp b/src/global/.global.cpp.swp deleted file mode 100644 index 80b64b8ba50f1d0cb3d76c0c5a98ae06f00f030b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 28672 zcmeI432-D=d4LBD1QW0^hT=4SmhJUOmNZ(M)81WME6uKjby?EdUUr?PHPbWG-k$DB zcaNmeVkjV$RJa`CB8Cu?0z0Hq6i1~vQ=!7)a0CK5On^WrAQYr9l}d3Ygh0Okb1s1AySZl5qN;S`M%tGLo3RTBjH=X#;WU0~EUC!|C zPy<5^oQ(#akh}V-{kp4PI&z7+@JYwcM#wPEPy<5^3^g#+z)%B24Gc9f)WA>!Lk;{N z)j-&IbnfNk{gHA0Ulsp-#E!p@ihtiUAl#$l>qh+df%tD*&ST==rTFjF@s0NVbK>jQ z$A5os$Nl*7@Xt^KLk$cyFx0?M149iAH89k`Py<5^3^g#+z)%B24V;Ar%yKSw5%H$Q z0U-PTv5W`L&*eS_ABDHW+u-%^Qn29|Ov9Bh2G4|Nz}L>peVt6+E^zpgeU&1@!K6oR%5)Q-la2;F%Km65P?w{bR@Co=MXhHx! zti$E-6u1x`1wVcq>BHCH%kU5IMc9BF;YBb3KRcJ}@BsKQ4F};Va3MSb{`9fA+&keO zxD#FhFNc>w3G(n{coIDLE4kdq;bU+Qgm61lUVv7 z4uftNJKNQzmRq(fIi(t=pG*Fz2NQ<^&Vg8l*9;rYFfb~%W1X0~O&wHMoSeGiTH+K6 zYRHjx$AX4-#vs@szfu+KlCHan`1$ zEvr(!JVAW95qesrqCYf_CBa)q%= z)X9pg&N7dp%6dR4qe4paipprL>ST-stG;geL0%WU-CXwRl8?t(=`Xa%TlYLiEjeWb zwP@AEZ>uj~^4uV#=0%6nfX7wx%7*DQximn6j$E+G<%L&bq@sYCrVw6L^yPd1^+x(p0;=e}C=6~9H6grOOZu=5`Oa8xM zivO1CR@F2XI^X6y-&Q)`Zsgmh1m2W{wQ_x{9MtTJBwlpOM)#*032EboW4YC^MuR4A z4c{vnWU2GLmH56rqVBj>hJ?i>+WT6r^KG$s>zed{{Mcy9#ggTR6PE8k_XQK278}jw zq<==eQw7pL(R-=U%F}RzsE-!L)Qs0?(FR>RsL=}vA^nubNARe72m7zQQX+5GWKhcI z?I38FB`ZIMj8V03pUMaO)d3Y;#V_f0$XxU7x^k>ar~cuHs zH8GK%g2>bo)452V6fGVif}W0k)CvPKQC40>g+61N9`BIi zFapST^cPxMbv>P1i|v&f)tl#qN|EO&$)hXnMSksR8auQa4#eDJ0iSc+3}imY}J&r%49P<~pN<);SYcbYfB=e8RU7 zO@s~8x6O6O3Rs{qSo>)8TGUz?;}qv9i6>XILaiVm;}S?Kn^1)@w<#l}1QV4L!qQV? zoDQ0+$5Emhq0f%R0hP$-rgT%8TN6cUXK6>-bdTPx1S*f97W0jTI!X?I%c^GaYuzfD z&1epfZg*8wR+h{$Qqj?Hkesrk#&n0W93r2lU3Er_3Ej+OQr)6!Pf4#ktyQyHb@NPV zj@s&J{0k;ot5cV(a;VoTvs&ErCYb?cc@RXG(d1T7^A~+1v}dde?YiYph8sp;m1M-2 zoU$gKCoQUdox+x+A5^VSu6dQR*~;(Jx5uLefZUB|8M(~76UCM(zw9+>5Q+|KtZ^?d z(`!#@+BuebD>Z%JB9E-EN~g4pp6sW~(j7SNkJGEG#g(;u)vTN9kd%cX|EO8#SG~zh zNq4NW?qQgb4&y4%wI-zlQaFimQ2x|S!`AC7+LE+(eTjnTb-^md(DH()V%p(*N<^S2 z7X5)4n32Kr)rw7(xeN}TvqADS)}95&V8PwRZ9Mo$B`zAX(`(bla&g62o?cm-n_#!GFt42*DQz*2fwHM=k5%pp!Uj21NP{=4| z*D5eDBy!#82WX0idVaTpf^@dY+M#49y`tE8jgqbRU6UdtqcWh>dasniRzu70CLn5y zM$Ps~`PvAJm+lHY)_QwDii*V)BGZm;^Q_giofApX2kpiFKL=a#nb?|Q|6d>5dY{9V zzZYHuCtw^dg=fKekb}>km&^SzycO<-*TObzKou+~!w6gq-@>;4I=lhyghTLj_z|}J z7vLnEfE(b6a4tL+K8fAF4J{ah?_-}o2ycbk;c4(x_ysokU&9Hw1ulTc!H==OKMgO2 zo8f720sI$b{#z)+OW-oN2t*gpfZ=DTfuRP58W?Ke|D^`*=&t@)PPJW|rf*ank44%+ z?Hruc%S^-b^CRN;Y~RJjF*P}*SGSQJ(q5G0M{b#3S)5xudQEq^5jn=no)x%PgskkM zRZ>!1H|(UEwJK)Q(N0D@i8^^%giT$72T>;|#*A2!;trjxnSnv>&H9e~Y5zjJ{9B!w zo}XSZj?6DDt@M?-Y`K9Qwgwft=k~);>K=IyOR;+-jwH)HTU=b7TRT3W;32MP%QwYG zIjHb`K^_JL^vQf!l+Y)&M-gkq1+1CrwPP#ASQ%Igy;fb^1GLT6vF@T@x_Qa6`b>6S zkBozu%chgpjZ`<{1CfP!C_kdbl0mTzS)I2%^i8RXaaA{4n39`T+tq^`VncjyrmT0R zXLS3-a|>5%drAuwKVO^!*9bzZ5yYDm{XMz&zErBtOB|c64T;`1Ka%4>w+wM0dZaFc zx$gN`Z*YE$=3QF9JwbLF-)}}FtjiW2mH}{uNwk2d3b+>c@N3}WJx!6j(xE)nHl^Uqk8)z#v&{hlXc zjoei|EPK)tH(A&;9j7H;rH&()5S~+>2xOn5*W=|kE$P=i-noj2Urt|bd+0i?vGb+Vd1;ZhYxuaDo3b%t zhMt|*14uO6MiOKWls$lYBOmR}j3lFt;F+g~MK1QZ<<1^kUBlMa8;;S0*qsU!nLwy3 zI}34B4o+E{_*}}iKOPM`<7#30*7@S%?wRREOUugAv9;x6YsTu0#aoIit0=2ecBf?; z%cqR8;doW!R8;!PUP|8^F_Ea){}-^ltIK(N?Em)p{tsc_zX>h_vHO37{r*FE1Bh*3 zf#<{X;B(mUe+8d`JK!34KKvWD`}g5@;del6`xn6Z@NMk)`{8boJ^Lak_*d-sPr>Wp z*I+-K4`0KMe;=HN1rXc*8`$voK>*K($H3>Y+usIv!$}y2XTyW|I6eU%f|tULa6WvU zee<`&ufsYN;K}e7_RsV16ZXhI0v;^EQ{m_Ai~kAy2HXUr@Jsaf4fp`u4Wds6dOkG! zLCtf!ZnJz#G|Nt%x44s`rX_f~psj3zBzU@@_VgHSN9iuOr$*QlFI}|5q*ju}rWwlY zOuwS&VZ9-P!_Ema&-GYgjH>Nl`$w59QpB`<(EAmErzWG(Td}wZMSg~Xq+fP=QOG_g zAPm#lG+Q!IWiLgqgaHL$0*dxdJ3FR0*~{iS+j3!Rf}N*^V^>-evXRZ+)P%Lw@Nt;) zA>v%s^?Yl=Ecy6x^^S9UQN;?2t%ow|VX0%5C)+W~O%70aIMO+=v1!?UndO+oMl3FD zy#W{EhMr``VNa}7v)D?&1i|^w-b~pJWP>W+%MnMn>uoZHh|R;cgXvZ$Lc82j`(Gf_ zX*Ov-1}MlxlcVV!XHcF|K?i6#+E z9$McEH_U3`@#H}au`JB+>Z|ZLS6vH#zai_Ez0PUr0upH}KD#9*>z{{|7)zzJF)KA= zvW|L4iLvDBBo@!q{Tl7*a~YEB#!IH64?5|(n=KBix}#=LbS=Mnw#w|pPbM5KL`#U6;WZ2ocV^kiI_ zCEJ5Yh6H-z4H6wVcFD>lRL2-3&i6t&s2~S|ZnvV-8halQSn>BoAC%BpCeV#GP;k5@ zn@J~(zj_~ckkDGI%mkW66FSBKae99*gD?rdL(~ovM|PUijGKNFOG5U@VkcN1y3!MK z$zy+;!-b3}eoJP26w^M086RL^lk)Wc+>yJEHGDEUDJFa4EN3W(tP?u9QJZoSm9vT3 zM3z<1BqulJ;D|>!f71>umP2Fw>l}St?z4-<#e*HKg{pN{E}WL)JP_s8*0Z)(m9|`?!&o}PE|xyW3B{roZ&@fFGC6HQ*5w=$YohpMocP*A@lxE1_NK9r;}Z<@YqG+j|r|h0qnkf>8VwWTxw~7}r^p+pe`KaYWG_cLir3#Ot7|BRT-e;oDQu zk>3h$oAWkf&zsJC!wdpLTBB_RM(A;9QQICO)QWxP6*tPf*(4X7Pd%tk?7M9*CqFrJ zcSz~Oqo>>2Pd`N(c_-rPC22Xv%8|&ZE(hvN*5f3_ro%|1!K6?t)XW2>ak%_!f5mCqef4UjfJA06ZCfiv9mLAU=WL z1{%oGfa1s0hU%&(Kes~QyFb$W$Pw@wQ8Qu+YzTSc(a1s0xpTHO4 z1Mpg?!wg&o|AC+2EAS3@85H5E@O69z?}bg6gQvj{@D+R%-V9Z^9?paBqW?b!sfV{n zUGPzSbd#7k)UZSRYO zFNZSLGV#(ymu<6=rSs_U>7ovA8{(0|uG~}RB=`V9t!#n_K0we`HbDeW7qpy74x!V< zY-H1r1Ro$Mn}#+9YG@;yh9r1;j%Z++bfm+li=tsT+1imaKu{~2Ac7ANw3ST|!P5mT zXOctcbTKrHEEsNZFU&@kbcjFTAb|et5$G4k?2#xf* zr~T(R#sG0nW^qxBbaAcBxr<{A5SKZ3t=)3Bl{t5DjP%@XWzJn3BVF7!&&*}$E)qOl z5a%ee2$JCGg0`{=lHlorWRH1g(Zv48oX)*lY){xXH#Ff0jKB}D`Q_aId*H25fdc#-d;bA=JDi49cpjV!-^JelGk6o6f;o5= z{0v+F@8JFLI%vQNxE3A<-^SkmDEtMy6K;nC@N{?(oBz}B`>+WOcp*Fi{vF%@es}}C z8UnZu9uJR(AL9dfAA~Ri7s6Mm6RDSX#r}XEkAR%$h>wwQnvKUXvv}(JBkyQ2JBDYI zN0D2Bo%x(gBHrE-F#6vfkZIsxLPK#c8klS%Yuz_^Bi;Ul0k|C4rx1FRmr1>nSJ8w%# z7m~&5=4~sfLU;j71}(|!T2jR@Rq%lL-p5HDK3$Z4V@w87(W_=sbH#gPveX!Fok?bnVj8#1EzTN8 zipAN(EbB;nmrS@kO2+Yv?72up9F&dK#cn2IRdq6zn3%gy{wKP5!*{E^k%BV|y4~$* z=;*{rWGGHIHH}Skqy3to3{5f-E0Inj>&wQvWAd_zy?2;6M!LA Date: Mon, 24 Jun 2024 13:03:29 -0700 Subject: [PATCH 40/62] Comment format for Cosmology::dtda_cosmo() Co-authored-by: Matthew Abruzzo --- src/cosmology/cosmology_functions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 7a1ee8afe..02d48d367 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -24,9 +24,9 @@ void Grid3D::Initialize_Cosmology(struct Parameters *P) chprintf("Cosmology Successfully Initialized. \n\n"); } +/* Computes dt/da * da */ Real Cosmology::dtda_cosmo(Real da, Real a) { - // Return dt/da * da Real a2 = a * a; Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - current_a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / a + a2 * Omega_L * fac_de + Omega_K) * H0; From 6063e44dafe66e4da2bacef303fe1f5ff0870420 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 13:04:04 -0700 Subject: [PATCH 41/62] Comment format for Cosmology::Get_dt_from_da_rk Co-authored-by: Matthew Abruzzo --- src/cosmology/cosmology_functions.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 02d48d367..573bdd605 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -33,10 +33,11 @@ Real Cosmology::dtda_cosmo(Real da, Real a) return da / a_dot; } + +/* Compute dt/da * da. dt/da is computed with a Runge-Kutta integration step */ Real Cosmology::Get_dt_from_da_rk(Real da, Real a) { - // Return dt/da * da - // But compute dt/da using a Runge-Kutta integration step + Real a3 = 0.3; Real a4 = 0.6; Real a5 = 1.0; From 8df15e0561acd8251daec6bb456a645bf8caf5d2 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 13:06:40 -0700 Subject: [PATCH 42/62] Comment format for Cosmology::Create_Expansion_History_File() Co-authored-by: Matthew Abruzzo --- src/cosmology/cosmology_functions.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 573bdd605..7d67b14a0 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -186,9 +186,9 @@ static inline bool Is_Root_Proc() #endif } +/* create the file for recording the expansion history */ void Cosmology::Create_Expansion_History_File(struct Parameters *P) { - // writing the expansion history if (not Is_Root_Proc()) { return; } From 53beabb999a554d424b872f2c2b1a2bfb51c90f3 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 13:07:40 -0700 Subject: [PATCH 43/62] Comment format for Cosmology::Write_Expansion_History_Entry() Co-authored-by: Matthew Abruzzo --- src/cosmology/cosmology_functions.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index 7d67b14a0..ea4ab6a27 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -221,6 +221,7 @@ void Cosmology::Create_Expansion_History_File(struct Parameters *P) out_file.close(); } +/* Write the current entry to the expansion history file */ void Cosmology::Write_Expansion_History_Entry(void) { if (not Is_Root_Proc()) { From 2b3b2a0f08be0754c7e747f447a0f23ae3ff485b Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 16:26:48 -0700 Subject: [PATCH 44/62] Adding requested comments Comments about the universal time offset and temporary retained but commented code. --- src/cosmology/cosmology.cpp | 15 +++++++++++++-- src/cosmology/cosmology_functions.cpp | 6 ++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index 06a69f794..25de5a4a3 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -48,16 +48,27 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part time_conversion = KPC; t_secs = 0; - // get time offset + // The following code computes the universal time + // at the scale factor of the ICs or restart file + // Pick a small scale factor step for integrating the universal time Real da_t_sec = 1.0e-2 * current_a; - // small but non-zero + // Pick a small but non-zero starting scale factor for the integral Real a_t_sec = 1.0e-6; + // Step for the time integral, corresponding to da_t_sec Real dt_physical; + + // Advance a_t_sec until it matches current_a, and integrate time while (a_t_sec < current_a) { + + // Limit the scale a_t_sec factor to current_a if (a_t_sec + da_t_sec > current_a) { da_t_sec = current_a - a_t_sec; } + + // Compute the time step dt_physical = Get_dt_from_da(da_t_sec, a_t_sec); + + // Advance the time in seconds and the scale factor t_secs += dt_physical * time_conversion; a_t_sec += da_t_sec; // chprintf(" Revised a_t_sec %f da_t_sec %f t_start : %e Myr\n", a_t_sec, da_t_sec, t_secs / MYR); diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index ea4ab6a27..f2af5b308 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -71,6 +71,12 @@ Real Cosmology::Get_da_from_dt(Real dt) Real Cosmology::Get_dt_from_da(Real da, Real a) { return Get_dt_from_da_rk(da, a); + + /* The following commented code was the original Euler + integrator for computing time from the scale factor. + This has been left here temporarily to ease comparison + with the Runge-Kutta integrator, but it can be removed + eventually. */ /* Real a2 = a * a; Real fac_de = pow(a, -3 * (1 + w0 + wa)) * exp(-3 * wa * (1 - a)); Real a_dot = sqrt(Omega_R / a2 + Omega_M / a + a2 * Omega_L * fac_de + Omega_K) * H0; From 09c9fbb6ce83235809070b653e8dc5b6869dcee9 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 16:29:58 -0700 Subject: [PATCH 45/62] Changes to frontier make and setup to compile w/ Following the suggestion of Trey White, confirmed on Frontier --- builds/make.host.frontier | 4 ++-- builds/setup.frontier.cce.sh | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/builds/make.host.frontier b/builds/make.host.frontier index bae874c78..92be53577 100644 --- a/builds/make.host.frontier +++ b/builds/make.host.frontier @@ -5,9 +5,9 @@ CXX = CC GPUCXX ?= hipcc CXXFLAGS_DEBUG = -g -O0 -std=c++17 -CXXFLAGS_OPTIMIZE = -g -Ofast -std=c++17 -Wno-unused-result +CXXFLAGS_OPTIMIZE = -g -Ofast -std=c++17 -Wno-unused-result --gcc-toolchain=${GCC_PATH}/snos -GPUFLAGS_OPTIMIZE = -std=c++17 --offload-arch=gfx90a -Wall -Wno-unused-result +GPUFLAGS_OPTIMIZE = -std=c++17 --offload-arch=gfx90a -Wall -Wno-unused-result --gcc-toolchain=${GCC_PATH}/snos GPUFLAGS_DEBUG = -g -O0 -std=c++17 --offload-arch=gfx90a -Wall -Wno-unused-result HIPCONFIG = -I$(ROCM_PATH)/include $(shell hipconfig -C) # workaround for Rocm 5.2 warnings #HIPCONFIG = $(shell hipconfig -C) diff --git a/builds/setup.frontier.cce.sh b/builds/setup.frontier.cce.sh index 86999a87b..dfed5dbbd 100755 --- a/builds/setup.frontier.cce.sh +++ b/builds/setup.frontier.cce.sh @@ -9,6 +9,7 @@ module load rocm/5.7.1 module load craype-accel-amd-gfx90a module load cray-hdf5 cray-fftw module load googletest/1.10.0 +module load gcc-mixed #-- GPU-aware MPI export MPICH_GPU_SUPPORT_ENABLED=1 From 5d3c10bab9a0784e471aecb0fff893fbed548885 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 24 Jun 2024 16:47:29 -0700 Subject: [PATCH 46/62] Moves definition of Is_Root_Proc() to io.h Removes the duplicated function definition from cosmology_functions.cpp. --- src/cosmology/cosmology_functions.cpp | 11 ----------- src/io/io.cpp | 13 ++----------- src/io/io.h | 11 +++++++++++ 3 files changed, 13 insertions(+), 22 deletions(-) diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index f2af5b308..f735258ce 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -181,17 +181,6 @@ void Grid3D::Change_GAS_Frame_System(bool forward) } } -/* local function that designates whether we are using a root-process. It gives - * gives a sensible result regardless of whether we are using MPI */ -static inline bool Is_Root_Proc() -{ - #ifdef MPI_CHOLLA - return procID == root; - #else - return true; - #endif -} - /* create the file for recording the expansion history */ void Cosmology::Create_Expansion_History_File(struct Parameters *P) { diff --git a/src/io/io.cpp b/src/io/io.cpp index 12f18fe9c..87cc8a34c 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -36,17 +36,7 @@ * output routine */ void Rotate_Point(Real x, Real y, Real z, Real delta, Real phi, Real theta, Real *xp, Real *yp, Real *zp); -/* local function that designates whether we are using a root-process. It gives - * gives a sensible result regardless of whether we are using MPI */ -static inline bool Is_Root_Proc() -{ -#ifdef MPI_CHOLLA - return procID == root; -#else - return true; -#endif -} - +/* Generate the log output file */ void Create_Log_File(struct Parameters P) { if (not Is_Root_Proc()) { @@ -75,6 +65,7 @@ void Create_Log_File(struct Parameters P) out_file.close(); } +/* Write an entry in the log output file */ void Write_Message_To_Log_File(const char *message) { if (not Is_Root_Proc()) { diff --git a/src/io/io.h b/src/io/io.h index d8f6ca8ca..87ef86495 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -7,6 +7,17 @@ #include "../global/global.h" #include "../grid/grid3D.h" +/* Local function that designates whether we are using a root-process. It gives + * * gives a sensible result regardless of whether we are using MPI */ +static inline bool Is_Root_Proc() +{ +#ifdef MPI_CHOLLA + return procID == root; +#else + return true; +#endif +} + /* Write the data */ void Write_Data(Grid3D& G, struct Parameters P, int nfile); From 00d6f94f3b18b5b43e82f5c028fe15279a13b0cb Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Tue, 25 Jun 2024 20:44:31 -0700 Subject: [PATCH 47/62] Synching changes from compile-chemistry-gpu --- src/chemistry_gpu/chemistry_functions.cpp | 10 +-- src/chemistry_gpu/chemistry_functions_gpu.cu | 71 ++++++++++---------- src/chemistry_gpu/chemistry_gpu.h | 10 ++- 3 files changed, 50 insertions(+), 41 deletions(-) diff --git a/src/chemistry_gpu/chemistry_functions.cpp b/src/chemistry_gpu/chemistry_functions.cpp index 7999a6d55..3ce7dfaf6 100644 --- a/src/chemistry_gpu/chemistry_functions.cpp +++ b/src/chemistry_gpu/chemistry_functions.cpp @@ -59,7 +59,7 @@ void Grid3D::Initialize_Chemistry(struct Parameters *P) Chem.H.density_units = Chem.H.density_units / Chem.H.a_value / Chem.H.a_value / Chem.H.a_value; Chem.H.length_units = Chem.H.length_units / Cosmo.cosmo_h * Chem.H.a_value; Chem.H.time_units = Chem.H.time_units / Cosmo.cosmo_h; - Chem.H.dens_number_conv = Chem.H.density_number_conv * pow(Chem.H.a_value, 3); + Chem.H.dens_number_conv = Chem.H.dens_number_conv * pow(Chem.H.a_value, 3); #endif // COSMOLOGY Chem.H.velocity_units = Chem.H.length_units / Chem.H.time_units; @@ -74,10 +74,12 @@ void Grid3D::Initialize_Chemistry(struct Parameters *P) time_base = Chem.H.time_units; Chem.H.cooling_units = (pow(length_base, 2) * pow(MH, 2)) / (dens_base * pow(time_base, 3)); Chem.H.reaction_units = MH / (dens_base * time_base); - // printf(" cooling_units: %e\n", Chem.H.cooling_units ); - // printf(" reaction_units: %e\n", Chem.H.reaction_units ); + Chem.H.max_iter = 10000; - Chem.H.max_iter = 10000; + // The chemistry GPU functions need access to the temperature floor + // but they don't have access to P. Use Chem.H as a carrier. + // The P->temperature_floor is always defined, so safe to set here. + Chem.H.temperature_floor = P->temperature_floor; // Initialize all the rates Chem.Initialize(P); diff --git a/src/chemistry_gpu/chemistry_functions_gpu.cu b/src/chemistry_gpu/chemistry_functions_gpu.cu index 7c9bfe2cf..23fdff9f7 100644 --- a/src/chemistry_gpu/chemistry_functions_gpu.cu +++ b/src/chemistry_gpu/chemistry_functions_gpu.cu @@ -90,7 +90,7 @@ class Thermal_State } }; -__device__ void get_temperature_indx(Real T, Chemistry_Header &Chem_H, int &temp_indx, Real &delta_T, Real temp_old, +__device__ void get_temperature_indx(Real T, ChemistryHeader &Chem_H, int &temp_indx, Real &delta_T, Real temp_old, bool print) { Real logT, logT_start, d_logT, logT_l, logT_r; @@ -118,7 +118,7 @@ __device__ Real interpolate_rate(Real *rate_table, int indx, Real delta) return rate_val; } -__device__ Real Get_Cooling_Rates(Thermal_State &TS, Chemistry_Header &Chem_H, Real dens_number_conv, Real current_z, +__device__ Real Get_Cooling_Rates(Thermal_State &TS, ChemistryHeader &Chem_H, Real dens_number_conv, Real current_z, Real temp_prev, float photo_h_HI, float photo_h_HeI, float photo_h_HeII, bool print) { int temp_indx; @@ -201,7 +201,7 @@ __device__ Real Get_Cooling_Rates(Thermal_State &TS, Chemistry_Header &Chem_H, R return U_dot; } -__device__ void Get_Reaction_Rates(Thermal_State &TS, Chemistry_Header &Chem_H, Real &k_coll_i_HI, Real &k_coll_i_HeI, +__device__ void Get_Reaction_Rates(Thermal_State &TS, ChemistryHeader &Chem_H, Real &k_coll_i_HI, Real &k_coll_i_HeI, Real &k_coll_i_HeII, Real &k_coll_i_HI_HI, Real &k_coll_i_HI_HeI, Real &k_recomb_HII, Real &k_recomb_HeII, Real &k_recomb_HeIII, bool print) { @@ -257,7 +257,7 @@ __device__ Real linear_interpolation(Real delta_x, int indx_l, int indx_r, float return v; } -__device__ void Get_Current_UVB_Rates(Real current_z, Chemistry_Header &Chem_H, float &photo_i_HI, float &photo_i_HeI, +__device__ void Get_Current_UVB_Rates(Real current_z, ChemistryHeader &Chem_H, float &photo_i_HI, float &photo_i_HeI, float &photo_i_HeII, float &photo_h_HI, float &photo_h_HeI, float &photo_h_HeII, bool print) { @@ -287,7 +287,7 @@ __device__ void Get_Current_UVB_Rates(Real current_z, Chemistry_Header &Chem_H, photo_h_HeII = linear_interpolation(delta_x, indx_l, indx_l + 1, Chem_H.photo_heat_HeII_rate_d); } -__device__ Real Get_Chemistry_dt(Thermal_State &TS, Chemistry_Header &Chem_H, Real &HI_dot, Real &e_dot, Real U_dot, +__device__ Real Get_Chemistry_dt(Thermal_State &TS, ChemistryHeader &Chem_H, Real &HI_dot, Real &e_dot, Real U_dot, Real k_coll_i_HI, Real k_coll_i_HeI, Real k_coll_i_HeII, Real k_coll_i_HI_HI, Real k_coll_i_HI_HeI, Real k_recomb_HII, Real k_recomb_HeII, Real k_recomb_HeIII, float photo_i_HI, float photo_i_HeI, float photo_i_HeII, int n_iter, Real HI_dot_prev, @@ -328,7 +328,9 @@ __device__ Real Get_Chemistry_dt(Thermal_State &TS, Chemistry_Header &Chem_H, Re } #ifdef TEMPERATURE_FLOOR - if (TS.get_temperature(Chem_H.gamma) < TEMP_FLOOR) TS.U = TS.compute_U(TEMP_FLOOR, Chem_H.gamma); + if (TS.get_temperature(Chem_H.gamma) < Chem_H.temperature_floor) { + TS.U = TS.compute_U(Chem_H.temperature_floor, Chem_H.gamma); + } #endif energy = fmax(TS.U * TS.d, tiny); @@ -355,7 +357,7 @@ __device__ Real Get_Chemistry_dt(Thermal_State &TS, Chemistry_Header &Chem_H, Re return dt; } -__device__ void Update_Step(Thermal_State &TS, Chemistry_Header &Chem_H, Real dt, Real U_dot, Real k_coll_i_HI, +__device__ void Update_Step(Thermal_State &TS, ChemistryHeader &Chem_H, Real dt, Real U_dot, Real k_coll_i_HI, Real k_coll_i_HeI, Real k_coll_i_HeII, Real k_coll_i_HI_HI, Real k_coll_i_HI_HeI, Real k_recomb_HII, Real k_recomb_HeII, Real k_recomb_HeIII, float photo_i_HI, float photo_i_HeI, float photo_i_HeII, Real &HI_dot_prev, Real &e_dot_prev, Real &temp_prev, @@ -421,13 +423,15 @@ __device__ void Update_Step(Thermal_State &TS, Chemistry_Header &Chem_H, Real dt // Update internal energy TS.U += U_dot / TS.d * dt; #ifdef TEMPERATURE_FLOOR - if (TS.get_temperature(Chem_H.gamma) < TEMP_FLOOR) TS.U = TS.compute_U(TEMP_FLOOR, Chem_H.gamma); + if (TS.get_temperature(Chem_H.gamma) < Chem_H.temperature_floor) { + TS.U = TS.compute_U(Chem_H.temperature_floor, Chem_H.gamma); + } #endif if (print) printf("Updated U: %e \n", TS.U); } __global__ void Update_Chemistry_kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, - Real dt_hydro, Chemistry_Header Chem_H) + Real dt_hydro, ChemistryHeader Chem_H) { int id, xid, yid, zid, n_cells, n_iter; Real d, d_inv, vx, vy, vz; @@ -482,17 +486,16 @@ __global__ void Update_Chemistry_kernel(Real *dev_conserved, int nx, int ny, int dt_hydro = dt_hydro / Chem_H.time_units; #ifdef COSMOLOGY - dt_hydro *= current_a * current_a / Chem_H.H0 * 1000 * - KPC + dt_hydro *= current_a * current_a / Chem_H.H0 * 1000 * KPC; #endif // COSMOLOGY - // dt_hydro = dt_hydro * current_a * current_a / Chem_H.H0 * - // 1000 * KPC / Chem_H.time_units; - // delta_a = Chem_H.H0 * sqrt( Chem_H.Omega_M/current_a + - // Chem_H.Omega_L*pow(current_a, 2) ) / ( 1000 * KPC ) * - // dt_hydro * Chem_H.time_units; - - // Initialize the thermal state - Thermal_State TS; + // dt_hydro = dt_hydro * current_a * current_a / Chem_H.H0 * + // 1000 * KPC / Chem_H.time_units; + // delta_a = Chem_H.H0 * sqrt( Chem_H.Omega_M/current_a + + // Chem_H.Omega_L*pow(current_a, 2) ) / ( 1000 * KPC ) * + // dt_hydro * Chem_H.time_units; + + // Initialize the thermal state + Thermal_State TS; TS.d = dev_conserved[id] / a3; TS.d_HI = dev_conserved[id + n_cells * grid_enum::HI_density] / a3; TS.d_HII = dev_conserved[id + n_cells * grid_enum::HII_density] / a3; @@ -603,7 +606,7 @@ __global__ void Update_Chemistry_kernel(Real *dev_conserved, int nx, int ny, int } void Do_Chemistry_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, - Chemistry_Header &Chem_H) + ChemistryHeader &Chem_H) { float time; cudaEvent_t start, stop; @@ -648,7 +651,7 @@ void Do_Chemistry_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghos // Calculation of k1 (HI + e --> HII + 2e) // k1_rate -__device__ Real coll_i_HI_rate(Real T, Real units) +__host__ __device__ Real coll_i_HI_rate(Real T, Real units) { Real T_ev = T / 11605.0; Real logT_ev = log(T_ev); @@ -666,7 +669,7 @@ __device__ Real coll_i_HI_rate(Real T, Real units) // Calculation of k3 (HeI + e --> HeII + 2e) // k3_rate -__device__ Real coll_i_HeI_rate(Real T, Real units) +__host__ __device__ Real coll_i_HeI_rate(Real T, Real units) { Real T_ev = T / 11605.0; Real logT_ev = log(T_ev); @@ -684,7 +687,7 @@ __device__ Real coll_i_HeI_rate(Real T, Real units) // Calculation of k4 (HeII + e --> HeI + photon) // k4_rate -__device__ Real recomb_HeII_rate(Real T, Real units, bool use_case_B) +__host__ __device__ Real recomb_HeII_rate(Real T, Real units, bool use_case_B) { Real T_ev = T / 11605.0; Real logT_ev = log(T_ev); @@ -703,7 +706,7 @@ __device__ Real recomb_HeII_rate(Real T, Real units, bool use_case_B) } } // k4_rate Case A -__device__ Real recomb_HeII_rate_case_A(Real T, Real units) +__host__ __device__ Real recomb_HeII_rate_case_A(Real T, Real units) { Real T_ev = T / 11605.0; Real logT_ev = log(T_ev); @@ -716,7 +719,7 @@ __device__ Real recomb_HeII_rate_case_A(Real T, Real units) } } // k4_rate Case B -__device__ Real recomb_HeII_rate_case_B(Real T, Real units) +__host__ __device__ Real recomb_HeII_rate_case_B(Real T, Real units) { // If case B recombination on. return 1.26e-14 * pow(5.7067e5 / T, 0.75) / units; @@ -724,7 +727,7 @@ __device__ Real recomb_HeII_rate_case_B(Real T, Real units) // Calculation of k2 (HII + e --> HI + photon) // k2_rate -__device__ Real recomb_HII_rate(Real T, Real units, bool use_case_B) +__host__ __device__ Real recomb_HII_rate(Real T, Real units, bool use_case_B) { if (use_case_B) { if (T < 1.0e9) { @@ -750,7 +753,7 @@ __device__ Real recomb_HII_rate(Real T, Real units, bool use_case_B) } } // k2_rate Case A -__device__ Real recomb_HII_rate_case_A(Real T, Real units) +__host__ __device__ Real recomb_HII_rate_case_A(Real T, Real units) { if (T > 5500) { // Convert temperature to appropriate form. @@ -769,7 +772,7 @@ __device__ Real recomb_HII_rate_case_A(Real T, Real units) } // k2_rate Case B -__device__ Real recomb_HII_rate_case_B(Real T, Real units) +__host__ __device__ Real recomb_HII_rate_case_B(Real T, Real units) { if (T < 1.0e9) { return 4.881357e-6 * pow(T, -1.5) * pow((1.0 + 1.14813e2 * pow(T, -0.407)), -2.242) / units; @@ -780,7 +783,7 @@ __device__ Real recomb_HII_rate_case_B(Real T, Real units) // Calculation of k5 (HeII + e --> HeIII + 2e) // k5_rate -__device__ Real coll_i_HeII_rate(Real T, Real units) +__host__ __device__ Real coll_i_HeII_rate(Real T, Real units) { Real T_ev = T / 11605.0; Real logT_ev = log(T_ev); @@ -800,7 +803,7 @@ __device__ Real coll_i_HeII_rate(Real T, Real units) // Calculation of k6 (HeIII + e --> HeII + photon) // k6_rate -__device__ Real recomb_HeIII_rate(Real T, Real units, bool use_case_B) +__host__ __device__ Real recomb_HeIII_rate(Real T, Real units, bool use_case_B) { Real k6; // Has case B recombination setting. @@ -816,7 +819,7 @@ __device__ Real recomb_HeIII_rate(Real T, Real units, bool use_case_B) return k6; } // k6_rate Case A -__device__ Real recomb_HeIII_rate_case_A(Real T, Real units) +__host__ __device__ Real recomb_HeIII_rate_case_A(Real T, Real units) { Real k6; // Has case B recombination setting. @@ -824,7 +827,7 @@ __device__ Real recomb_HeIII_rate_case_A(Real T, Real units) return k6; } // k6_rate Case B -__device__ Real recomb_HeIII_rate_case_B(Real T, Real units) +__host__ __device__ Real recomb_HeIII_rate_case_B(Real T, Real units) { Real k6; // Has case B recombination setting. @@ -838,7 +841,7 @@ __device__ Real recomb_HeIII_rate_case_B(Real T, Real units) // Calculation of k57 (HI + HI --> HII + HI + e) // k57_rate -__device__ Real coll_i_HI_HI_rate(Real T, Real units) +__host__ __device__ Real coll_i_HI_HI_rate(Real T, Real units) { // These rate coefficients are from Lenzuni, Chernoff & Salpeter (1991). // k57 value based on experimental cross-sections from Gealy & van Zyl (1987). @@ -851,7 +854,7 @@ __device__ Real coll_i_HI_HI_rate(Real T, Real units) // Calculation of k58 (HI + HeI --> HII + HeI + e) // k58_rate -__device__ Real coll_i_HI_HeI_rate(Real T, Real units) +__host__ __device__ Real coll_i_HI_HeI_rate(Real T, Real units) { // These rate coefficients are from Lenzuni, Chernoff & Salpeter (1991). // k58 value based on cross-sections from van Zyl, Le & Amme (1981). diff --git a/src/chemistry_gpu/chemistry_gpu.h b/src/chemistry_gpu/chemistry_gpu.h index 79674c3a0..d9bbe0552 100644 --- a/src/chemistry_gpu/chemistry_gpu.h +++ b/src/chemistry_gpu/chemistry_gpu.h @@ -76,6 +76,10 @@ struct ChemistryHeader { float *photo_heat_HI_rate_d; float *photo_heat_HeI_rate_d; float *photo_heat_HeII_rate_d; + + // inherit the temperature floor + // from Parameters + float temperature_floor; }; #ifdef CHEMISTRY_GPU @@ -111,7 +115,7 @@ class Chem_GPU float *Ion_rates_HeI_d; float *Ion_rates_HeII_d; - struct Chemistry_Header H; + struct ChemistryHeader H; struct Fields { Real *temperature_h; @@ -152,7 +156,7 @@ n_ghost, int n_fields, Real dt, Real gamma) ionization fractions of H and He and update the internal energy to account for radiative cooling and photoheating from the UV background. */ void Do_Chemistry_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, - Chemistry_Header &Chem_H); + ChemistryHeader &Chem_H); #endif -#endif \ No newline at end of file +#endif From 196c8280593b9bab6a7908d19735f76db5366ec1 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Thu, 27 Jun 2024 21:08:29 -0700 Subject: [PATCH 48/62] Initialize cosmological chem+cooling ICs, don't read them The cosmological initial conditions by default do not have the chemistry fields recorded. While this change could be made, effectively at early times the universe is very neutral and at least for the current set of calculations we will initialize them directly. These changes allow for the cosmo chem+cooling runs to at least start. --- src/io/io.cpp | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 87cc8a34c..6b28644fe 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2304,6 +2304,7 @@ void Grid3D::Read_Grid(struct Parameters P) // close the file status = H5Fclose(file_id); #endif // BINARY or HDF5 + } /*! \fn Read_Grid_Binary(FILE *fp) @@ -2541,17 +2542,39 @@ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.dust_density, "/dust_density"); #endif // DUST - #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) + #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) + #ifdef COSMOLOGY + if (P.nfile > 0) { + #endif // COSMOLOGY Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HI_density, "/HI_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HII_density, "/HII_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeI_density, "/HeI_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeII_density, "/HeII_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeIII_density, "/HeIII_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.e_density, "/e_density"); - #ifdef GRACKLE_METALS + #ifdef GRACKLE_METALS Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.metal_density, "/metal_density"); - #endif // GRACKLE_METALS - #endif // COOLING_GRACKLE , CHEMISTRY_GPU + #endif // GRACKLE_METALS + #ifdef COSMOLOGY + } else { + // Initialize the density field + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; + C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; + C.HeI_density[id] = INITIAL_FRACTION_HEI * C.density[id]; + C.HeII_density[id] = INITIAL_FRACTION_HEII * C.density[id]; + C.HeIII_density[id] = INITIAL_FRACTION_HEIII * C.density[id]; + C.e_density[id] = INITIAL_FRACTION_ELECTRON * C.density[id]; + } + } + } + } + #endif // COSMOLOGY + #endif // COOLING_GRACKLE , CHEMISTRY_GPU #endif // SCALAR From 200d968447d156e5abddcfb5a29381e2ab888d9f Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 28 Jun 2024 21:36:54 -0700 Subject: [PATCH 49/62] Adding script to print initial statistics. This checks the cosmological simulation grid units. Something similar is in Bruno's original cosmology version. --- src/grid/grid3D.h | 6 ++ src/io/io.cpp | 214 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 220 insertions(+) diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index 4a1fa17a5..c77a6ebe2 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -530,6 +530,12 @@ class Grid3D void Read_Grid_HDF5(hid_t file_id, struct Parameters P); #endif +#if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) + /*! \fn void Print_Grid_Stats(struct Parameters P) + * \brief Compute stats for a grid property. */ + void Print_Grid_Stats(struct Parameters P); +#endif + /*! \fn void Reset(void) * \brief Reset the Grid3D class. */ void Reset(void); diff --git a/src/io/io.cpp b/src/io/io.cpp index 6b28644fe..ce05fcb52 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2497,6 +2497,220 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H dataset_buffer, grid_buffer); } +#if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) +/*! \fn void Print_Grid_Stats(struct Parameters P) + * \brief Compute stats for a grid property. */ +void Grid3D::Print_Grid_Stats(struct Parameters P) +{ + + // Print several interesting numbers + + // Density stats + mean_l = 0; + min_l = 1e65; + max_l = -1; + // Do density first + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += C.density[id]; + max_l = std::max(max_l, C.density[id]); + min_l = std::min(min_l, C.density[id]); + } + } + } + mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf( + " Density Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3] \n", + mean_l, min_l, max_l); + + // Momentum stats + + // x momenta + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += std::abs(momentum_x[id]); + max_l = std::max(max_l, std::abs(C.momentum_x[id])); + min_l = std::min(min_l, std::abs(C.momentum_x[id])); + } + } + } + mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf( + " abs(Momentum X) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", + mean_l, min_l, max_l); + + + // y momenta + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += std::abs(C.momentum_y[id]); + max_l = std::max(max_l, std::abs(C.momentum_y[id])); + min_l = std::min(min_l, std::abs(C.momentum_y[id])); + } + } + } + mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf( + " abs(Momentum Y) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", + mean_l, min_l, max_l); + + + // z momenta + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += std::abs(C.momentum_z[id]); + max_l = std::max(max_l, std::abs(C.momentum_z[id])); + min_l = std::min(min_l, std::abs(C.momentum_z[id])); + } + } + } + mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf( + " abs(Momentum Z) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", + mean_l, min_l, max_l); + + // Energy + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += C.Energy[id]; + max_l = std::max(max_l, C.Energy[id]); + min_l = std::min(min_l, C.Energy[id]); + } + } + } + mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf( + " Energy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ]\n", + mean_l, min_l, max_l); + + #ifdef DE + Real temp, temp_max_l, temp_min_l, temp_mean_l; + Real temp_min_g, temp_max_g, temp_mean_g; + temp_mean_l = 0; + temp_min_l = 1e65; + temp_max_l = -1; + mean_l = 0; + min_l = 1e65; + max_l = -1; + // Copy the internal Energy array to the grid + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real + 1; i++) { + id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + mean_l += C.GasEnergy[id]; + max_l = std::max(max_l, C.GasEnergy[id]); + min_l = std::min(min_l, C.GasEnergy[id]); + + temp = C.GasEnergy[id] / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; + temp_mean_l += temp; + temp_max_l = std::max(temp_max_l, C.GasEnergy[id]); + temp_min_l = std::min(temp_min_l, C.GasEnergy[id]); + } + } + } + mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); + temp_mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); + + + #if MPI_CHOLLA + mean_g = ReduceRealAvg( mean_l ); + max_g = ReduceRealMax( max_l ); + min_g = ReduceRealMin( min_l ); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + temp_mean_g = ReduceRealAvg( temp_mean_l ); + temp_max_g = ReduceRealMax( temp_max_l ); + temp_min_g = ReduceRealMin( temp_min_l ); + temp_mean_l = temp_mean_g; + temp_max_l = temp_max_g; + temp_min_l = temp_min_g; + #endif //MPI_CHOLLA + + chprintf( " GasEnergy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", + mean_l, min_l, max_l ); + chprintf( " Temperature Mean: %f Min: %f Max: %f [ K ] \n", + temp_mean_l, temp_min_l, temp_max_l ); + #endif // DE +} +#endif // PRINT_INITIAL_STATS and COSMOLOGY + /*! \fn void Read_Grid_HDF5(hid_t file_id) * \brief Read in grid data from an hdf5 file. */ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) From 8a35d5e8bd4a3d434dda4970822643fe8c521500 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 28 Jun 2024 21:44:13 -0700 Subject: [PATCH 50/62] Missing declarations and some typos. --- src/io/io.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index ce05fcb52..e683abd02 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2502,6 +2502,9 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H * \brief Compute stats for a grid property. */ void Grid3D::Print_Grid_Stats(struct Parameters P) { + int i, j, k, id, buf_id; + Real mean_l, min_l, max_l; + Real mean_g, min_g, max_g; // Print several interesting numbers @@ -2546,7 +2549,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) for (i = 0; i < H.nx_real + 1; i++) { id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); - mean_l += std::abs(momentum_x[id]); + mean_l += std::abs(C.momentum_x[id]); max_l = std::max(max_l, std::abs(C.momentum_x[id])); min_l = std::min(min_l, std::abs(C.momentum_x[id])); } From 3af012925cbb9b7155c03d60ec96b47d74de88ab Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 28 Jun 2024 22:08:59 -0700 Subject: [PATCH 51/62] Indexing error in print grid stats. --- src/io/io.cpp | 44 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index e683abd02..14570c7e6 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2513,10 +2513,12 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) min_l = 1e65; max_l = -1; // Do density first + /*for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny;*/ for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.density[id]; max_l = std::max(max_l, C.density[id]); @@ -2524,7 +2526,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) } } } - mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); @@ -2546,8 +2548,8 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_x[id]); max_l = std::max(max_l, std::abs(C.momentum_x[id])); @@ -2555,7 +2557,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) } } } - mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); @@ -2576,8 +2578,8 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_y[id]); max_l = std::max(max_l, std::abs(C.momentum_y[id])); @@ -2585,7 +2587,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) } } } - mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); @@ -2606,8 +2608,8 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_z[id]); max_l = std::max(max_l, std::abs(C.momentum_z[id])); @@ -2615,7 +2617,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) } } } - mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); @@ -2635,8 +2637,8 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.Energy[id]; max_l = std::max(max_l, C.Energy[id]); @@ -2644,7 +2646,7 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) } } } - mean_l /= ((H.nz_real + 1) * (H.ny_real) * (H.nx_real)); + mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); @@ -2667,14 +2669,10 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) mean_l = 0; min_l = 1e65; max_l = -1; - // Copy the internal Energy array to the grid - mean_l = 0; - min_l = 1e65; - max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.GasEnergy[id]; max_l = std::max(max_l, C.GasEnergy[id]); @@ -2777,8 +2775,8 @@ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) // Initialize the density field for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real + 1; i++) { - id = (i + H.n_ghost - 1) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; From 6ae91bb01a8ea0ae72a37762de79952e7093594d Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Fri, 28 Jun 2024 22:14:17 -0700 Subject: [PATCH 52/62] Replacing GasEnergy with temp in print grid stats --- src/io/io.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 14570c7e6..736e64856 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2680,8 +2680,8 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) temp = C.GasEnergy[id] / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; temp_mean_l += temp; - temp_max_l = std::max(temp_max_l, C.GasEnergy[id]); - temp_min_l = std::min(temp_min_l, C.GasEnergy[id]); + temp_max_l = std::max(temp_max_l, temp); + temp_min_l = std::min(temp_min_l, temp); } } } From f2da5b6dbaeb4449accae2d1f1336c659cc6355b Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Sun, 30 Jun 2024 22:12:04 -0700 Subject: [PATCH 53/62] Backporting cooling+chemistry bug fix Removing the dens_number_conv a^3 scaling from chemistry_functions.cpp --- src/chemistry_gpu/chemistry_functions.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/chemistry_gpu/chemistry_functions.cpp b/src/chemistry_gpu/chemistry_functions.cpp index 3ce7dfaf6..d6074e21d 100644 --- a/src/chemistry_gpu/chemistry_functions.cpp +++ b/src/chemistry_gpu/chemistry_functions.cpp @@ -55,11 +55,10 @@ void Grid3D::Initialize_Chemistry(struct Parameters *P) Chem.H.time_units = kpc_km; Chem.H.dens_number_conv = Chem.H.density_units / MH; #ifdef COSMOLOGY - Chem.H.a_value = Cosmo.current_a; - Chem.H.density_units = Chem.H.density_units / Chem.H.a_value / Chem.H.a_value / Chem.H.a_value; - Chem.H.length_units = Chem.H.length_units / Cosmo.cosmo_h * Chem.H.a_value; - Chem.H.time_units = Chem.H.time_units / Cosmo.cosmo_h; - Chem.H.dens_number_conv = Chem.H.dens_number_conv * pow(Chem.H.a_value, 3); + Chem.H.a_value = Cosmo.current_a; + Chem.H.density_units = Chem.H.density_units / Chem.H.a_value / Chem.H.a_value / Chem.H.a_value; + Chem.H.length_units = Chem.H.length_units / Cosmo.cosmo_h * Chem.H.a_value; + Chem.H.time_units = Chem.H.time_units / Cosmo.cosmo_h; #endif // COSMOLOGY Chem.H.velocity_units = Chem.H.length_units / Chem.H.time_units; From 6f84e2e74935040d0c8cafc3a1f73c1f9e85b9ef Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Sun, 30 Jun 2024 22:20:58 -0700 Subject: [PATCH 54/62] Functions for printing grid statistics. --- src/grid/grid3D.h | 6 ++-- src/io/io.cpp | 76 ++++++++++++++++++++++++++++++++++++++++++----- src/io/io.h | 4 +++ 3 files changed, 75 insertions(+), 11 deletions(-) diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index c77a6ebe2..7d63385ec 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -531,9 +531,9 @@ class Grid3D #endif #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) - /*! \fn void Print_Grid_Stats(struct Parameters P) - * \brief Compute stats for a grid property. */ - void Print_Grid_Stats(struct Parameters P); + /*! \fn void Print_Grid_Stats(void) + * \brief Compute stats for grid properties. */ + void Print_Grid_Stats(void) #endif /*! \fn void Reset(void) diff --git a/src/io/io.cpp b/src/io/io.cpp index 736e64856..63cae809e 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2498,9 +2498,20 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H } #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) -/*! \fn void Print_Grid_Stats(struct Parameters P) - * \brief Compute stats for a grid property. */ -void Grid3D::Print_Grid_Stats(struct Parameters P) + +/*! \fn void Print_Stats(void) + * \brief Compute stats for a grid. */ +void Print_Stats(Grid3D &G) +{ + // Synchronize + cudaMemcpy(G.C.density, G.C.device, G.H.n_fields * G.H.n_cells * sizeof(Real), cudaMemcpyDeviceToHost); + // Write data + G.Print_Grid_Stats(); +} + +/*! \fn void Print_Grid_Stats(void) + * \brief Compute stats for grid properties. */ +void Grid3D::Print_Grid_Stats(void) { int i, j, k, id, buf_id; Real mean_l, min_l, max_l; @@ -2513,8 +2524,6 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) min_l = 1e65; max_l = -1; // Do density first - /*for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny;*/ for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { @@ -2660,9 +2669,60 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) " Energy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ]\n", mean_l, min_l, max_l); - #ifdef DE Real temp, temp_max_l, temp_min_l, temp_mean_l; Real temp_min_g, temp_max_g, temp_mean_g; + Real gase, vx, vy, vz; + temp_mean_l = 0; + temp_min_l = 1e65; + temp_max_l = -1; + mean_l = 0; + min_l = 1e65; + max_l = -1; + for (k = 0; k < H.nz_real; k++) { + for (j = 0; j < H.ny_real; j++) { + for (i = 0; i < H.nx_real; i++) { + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + vx = C.momentum_x[id] / C.density[id]; + vy = C.momentum_y[id] / C.density[id]; + vz = C.momentum_z[id] / C.density[id]; + gase = C.Energy[id] - 0.5 * C.density[id] * (vx * vx + vy * vy + vz * vz); + mean_l += gase; + max_l = std::max(max_l, gase); + min_l = std::min(min_l, gase); + + temp = gase / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; + temp_mean_l += temp; + temp_max_l = std::max(temp_max_l, temp); + temp_min_l = std::min(temp_min_l, temp); + } + } + } + mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); + temp_mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); + + + #if MPI_CHOLLA + mean_g = ReduceRealAvg( mean_l ); + max_g = ReduceRealMax( max_l ); + min_g = ReduceRealMin( min_l ); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + temp_mean_g = ReduceRealAvg( temp_mean_l ); + temp_max_g = ReduceRealMax( temp_max_l ); + temp_min_g = ReduceRealMin( temp_min_l ); + temp_mean_l = temp_mean_g; + temp_max_l = temp_max_g; + temp_min_l = temp_min_g; + #endif //MPI_CHOLLA + + chprintf( " GasEnergyCalc Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", + mean_l, min_l, max_l ); + chprintf( " TemperatureCalc Mean: %f Min: %f Max: %f [ K ] \n", + temp_mean_l, temp_min_l, temp_max_l ); + + #ifdef DE temp_mean_l = 0; temp_min_l = 1e65; temp_max_l = -1; @@ -2704,9 +2764,9 @@ void Grid3D::Print_Grid_Stats(struct Parameters P) temp_min_l = temp_min_g; #endif //MPI_CHOLLA - chprintf( " GasEnergy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", + chprintf( " GasEnergyDE Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", mean_l, min_l, max_l ); - chprintf( " Temperature Mean: %f Min: %f Max: %f [ K ] \n", + chprintf( " TemperatureDE Mean: %f Min: %f Max: %f [ K ] \n", temp_mean_l, temp_min_l, temp_max_l ); #endif // DE } diff --git a/src/io/io.h b/src/io/io.h index 87ef86495..8e3d5e66e 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -18,6 +18,10 @@ static inline bool Is_Root_Proc() #endif } +/*! \fn void Print_Stats(void) + * \brief Compute stats for a grid. */ +void Print_Stats(Grid3D &G); + /* Write the data */ void Write_Data(Grid3D& G, struct Parameters P, int nfile); From 58f1997a5b6a44519033e3771d08d916ae72faef Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Sun, 30 Jun 2024 22:30:52 -0700 Subject: [PATCH 55/62] Missing semicolon --- src/grid/grid3D.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index 7d63385ec..3f5c5663d 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -533,7 +533,7 @@ class Grid3D #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) /*! \fn void Print_Grid_Stats(void) * \brief Compute stats for grid properties. */ - void Print_Grid_Stats(void) + void Print_Grid_Stats(void); #endif /*! \fn void Reset(void) From 5466c906cc3e9317cb85ca648144e04d746ff1de Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Sun, 30 Jun 2024 23:56:54 -0700 Subject: [PATCH 56/62] Changes to DE to mimic original cosmo code Default dev branch does not match original cosmo code behavior. Removing the eta_1 logic from hydro_cuda.cu and then changing the eta_1 default value to 0.001 seems to roughly follow Bruno's cosmo verison. --- src/global/global.h | 3 ++- src/hydro/hydro_cuda.cu | 12 ++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/global/global.h b/src/global/global.h index fa7ee1b29..a2fd9178a 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -72,7 +72,8 @@ typedef double Real; // - In the future, we run tests and revisit the choice of DE_ETA_1 in // cosmological simulations #ifdef COSMOLOGY - #define DE_ETA_1 10.0 +// #define DE_ETA_1 10.0 + #define DE_ETA_1 0.001 #else #define DE_ETA_1 \ 0.001 // Ratio of U to E for which Internal Energy is used to compute the diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu index 1b223708b..bd67e94b3 100644 --- a/src/hydro/hydro_cuda.cu +++ b/src/hydro/hydro_cuda.cu @@ -944,7 +944,11 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) + #ifdef COSMOLOGY + if (U_total / Emax > eta_2) { + #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { + #endif U = U_total; } else { U = U_advected; @@ -1006,7 +1010,11 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) + #ifdef COSMOLOGY + if (U_total / Emax > eta_2) { + #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { + #endif U = U_total; } else { U = U_advected; @@ -1075,7 +1083,11 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) + #ifdef COSMOLOGY + if (U_total / Emax > eta_2) { + #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { + #endif U = U_total; } else { U = U_advected; From 7e0fbb28551f122623d6dcdad2cd4c789cec4437 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 1 Jul 2024 14:01:44 -0700 Subject: [PATCH 57/62] Roll back DE to dev version -- unnecessary change. --- src/global/global.h | 3 +-- src/hydro/hydro_cuda.cu | 12 ------------ 2 files changed, 1 insertion(+), 14 deletions(-) diff --git a/src/global/global.h b/src/global/global.h index a2fd9178a..fa7ee1b29 100644 --- a/src/global/global.h +++ b/src/global/global.h @@ -72,8 +72,7 @@ typedef double Real; // - In the future, we run tests and revisit the choice of DE_ETA_1 in // cosmological simulations #ifdef COSMOLOGY -// #define DE_ETA_1 10.0 - #define DE_ETA_1 0.001 + #define DE_ETA_1 10.0 #else #define DE_ETA_1 \ 0.001 // Ratio of U to E for which Internal Energy is used to compute the diff --git a/src/hydro/hydro_cuda.cu b/src/hydro/hydro_cuda.cu index bd67e94b3..1b223708b 100644 --- a/src/hydro/hydro_cuda.cu +++ b/src/hydro/hydro_cuda.cu @@ -944,11 +944,7 @@ __global__ void Select_Internal_Energy_1D(Real *dev_conserved, int nx, int n_gho // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) - #ifdef COSMOLOGY - if (U_total / Emax > eta_2) { - #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { - #endif U = U_total; } else { U = U_advected; @@ -1010,11 +1006,7 @@ __global__ void Select_Internal_Energy_2D(Real *dev_conserved, int nx, int ny, i // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) - #ifdef COSMOLOGY - if (U_total / Emax > eta_2) { - #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { - #endif U = U_total; } else { U = U_advected; @@ -1083,11 +1075,7 @@ __global__ void Select_Internal_Energy_3D(Real *dev_conserved, int nx, int ny, i // We only use the "advected" internal energy if both: // - the thermal energy divided by total energy is a small fraction (smaller than eta_1) // - AND we aren't masking shock heating (details controlled by Emax & eta_2) - #ifdef COSMOLOGY - if (U_total / Emax > eta_2) { - #else if ((U_total / E > eta_1) or (U_total / Emax > eta_2)) { - #endif U = U_total; } else { U = U_advected; From 1f05b8db9ee7ba7d217cf5509fc6ad533f9dff6f Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 1 Jul 2024 14:13:07 -0700 Subject: [PATCH 58/62] Change logic in PPMP so SIMPLE avoids the CTU ifdefs Changes were made to PPMP that affect reconstruction using the SIMPLE integrator, where "#ifdef CTU" was changed to "#ifndef VL". This affects the cosmological simulation hydro evolution, and behavior similar to Bruno's original code is recovered with this change reversed and SIMPLE integration used. --- src/reconstruction/ppmp_cuda.cu | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/reconstruction/ppmp_cuda.cu b/src/reconstruction/ppmp_cuda.cu index 2038f215a..59efac774 100644 --- a/src/reconstruction/ppmp_cuda.cu +++ b/src/reconstruction/ppmp_cuda.cu @@ -63,8 +63,8 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou // declare other variables Real del_q_imo, del_q_i, del_q_ipo; - #ifndef VL - // #ifdef CTU + //#ifndef VL + #ifdef CTU Real cs, cl, cr; // sound speed in cell i, and at left and right boundaries Real del_d, del_vx, del_vy, del_vz, del_p; // "slope" accross cell i Real d_6, vx_6, vy_6, vz_6, p_6; @@ -83,8 +83,8 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou #ifdef DE Real ge_i, ge_imo, ge_ipo, ge_imt, ge_ipt, ge_L, ge_R, E_kin, E, dge; - #ifndef VL - // #ifdef CTU + // #ifndef VL + #ifdef CTU Real del_ge, ge_6, geL_0, geR_0; #endif // CTU #endif // DE @@ -92,8 +92,8 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou #ifdef SCALAR Real scalar_i[NSCALARS], scalar_imo[NSCALARS], scalar_ipo[NSCALARS], scalar_imt[NSCALARS], scalar_ipt[NSCALARS]; Real scalar_L[NSCALARS], scalar_R[NSCALARS]; - #ifndef VL - // #ifdef CTU + // #ifndef VL + #ifdef CTU Real del_scalar[NSCALARS], scalar_6[NSCALARS], scalarL_0[NSCALARS], scalarR_0[NSCALARS]; #endif // CTU #endif // SCALAR @@ -471,8 +471,8 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou #endif // SCALAR #endif // FLATTENING - #ifndef VL - // #ifdef CTU + // #ifndef VL + #ifdef CTU // compute sound speed in cell i cs = sqrt(gamma * p_i / d_i); From 6ee56dafb318af3938c94a12b3761009b3ea834f Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Mon, 1 Jul 2024 19:26:23 -0700 Subject: [PATCH 59/62] Added minor changes to compile Lya analysis. --- builds/make.host.lux | 4 ++-- src/analysis/io_analysis.cpp | 3 +++ src/analysis/lya_statistics.cpp | 1 + 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/builds/make.host.lux b/builds/make.host.lux index 538cdddd5..861fdfcca 100644 --- a/builds/make.host.lux +++ b/builds/make.host.lux @@ -11,8 +11,8 @@ OMP_NUM_THREADS = 10 #-- Library CUDA_ROOT = /cm/shared/apps/cuda11.2/toolkit/current HDF5_ROOT = /cm/shared/apps/hdf5/1.10.6 -FFTW_ROOT = /home/brvillas/code/fftw-3.3.8 -PFFT_ROOT = /data/groups/comp-astro/bruno/code_mpi_local/pfft +FFTW_ROOT = /data/groups/comp-astro/cholla/code/fftw-3.3.10 +PFFT_ROOT = /data/groups/comp-astro/cholla/code/pfft GRACKLE_ROOT = /home/brvillas/code/grackle #Paris does not do GPU_MPI transfers diff --git a/src/analysis/io_analysis.cpp b/src/analysis/io_analysis.cpp index 341d83a61..ac7892fe6 100644 --- a/src/analysis/io_analysis.cpp +++ b/src/analysis/io_analysis.cpp @@ -9,6 +9,9 @@ // #define OUTPUT_SKEWERS_TRANSMITTED_FLUX +using std::string; +using std::ifstream; + #ifdef OUTPUT_SKEWERS void Grid3D::Output_Skewers_File(struct Parameters *P) { diff --git a/src/analysis/lya_statistics.cpp b/src/analysis/lya_statistics.cpp index 968011bae..50b69b0b5 100644 --- a/src/analysis/lya_statistics.cpp +++ b/src/analysis/lya_statistics.cpp @@ -13,6 +13,7 @@ // #define PRINT_ANALYSIS_LOG +using std::max; void AnalysisModule::Transfer_Skewers_Global_Axis(int axis) { bool am_I_root; From 116cd7a120b72757d157c74c1da5b56eb4642850 Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Tue, 2 Jul 2024 13:51:38 -0700 Subject: [PATCH 60/62] Clang formatting --- src/analysis/io_analysis.cpp | 4 +- src/cosmology/cosmology.cpp | 3 +- src/grid/grid3D.h | 2 +- src/io/io.cpp | 155 ++++++++++++++++---------------- src/io/io.h | 5 +- src/reconstruction/ppmp_cuda.cu | 2 +- 6 files changed, 83 insertions(+), 88 deletions(-) diff --git a/src/analysis/io_analysis.cpp b/src/analysis/io_analysis.cpp index ac7892fe6..88a94ccce 100644 --- a/src/analysis/io_analysis.cpp +++ b/src/analysis/io_analysis.cpp @@ -9,8 +9,8 @@ // #define OUTPUT_SKEWERS_TRANSMITTED_FLUX -using std::string; -using std::ifstream; + using std::string; + using std::ifstream; #ifdef OUTPUT_SKEWERS void Grid3D::Output_Skewers_File(struct Parameters *P) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index 25de5a4a3..a124303e5 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -58,7 +58,8 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part Real dt_physical; // Advance a_t_sec until it matches current_a, and integrate time - while (a_t_sec < current_a) { + while (a_t_sec < current_a) + { // Limit the scale a_t_sec factor to current_a if (a_t_sec + da_t_sec > current_a) { diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h index 3f5c5663d..7b6a7ae59 100644 --- a/src/grid/grid3D.h +++ b/src/grid/grid3D.h @@ -534,7 +534,7 @@ class Grid3D /*! \fn void Print_Grid_Stats(void) * \brief Compute stats for grid properties. */ void Print_Grid_Stats(void); -#endif +#endif /*! \fn void Reset(void) * \brief Reset the Grid3D class. */ diff --git a/src/io/io.cpp b/src/io/io.cpp index 63cae809e..236102fa8 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2304,7 +2304,6 @@ void Grid3D::Read_Grid(struct Parameters P) // close the file status = H5Fclose(file_id); #endif // BINARY or HDF5 - } /*! \fn Read_Grid_Binary(FILE *fp) @@ -2498,10 +2497,9 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H } #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) - -/*! \fn void Print_Stats(void) +/*! \fn void Print_Stats(Grid3D& G) * \brief Compute stats for a grid. */ -void Print_Stats(Grid3D &G) +void Print_Stats(Grid3D& G) { // Synchronize cudaMemcpy(G.C.density, G.C.device, G.H.n_fields * G.H.n_cells * sizeof(Real), cudaMemcpyDeviceToHost); @@ -2521,17 +2519,17 @@ void Grid3D::Print_Grid_Stats(void) // Density stats mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; // Do density first for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.density[id]; - max_l = std::max(max_l, C.density[id]); - min_l = std::min(min_l, C.density[id]); + max_l = std::max(max_l, C.density[id]); + min_l = std::min(min_l, C.density[id]); } } } @@ -2539,30 +2537,29 @@ void Grid3D::Print_Grid_Stats(void) #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; + max_l = max_g; + min_l = min_g; #endif // MPI_CHOLLA - chprintf( - " Density Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3] \n", - mean_l, min_l, max_l); + chprintf("Density Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3] \n", + mean_l, min_l, max_l); // Momentum stats // x momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_x[id]); - max_l = std::max(max_l, std::abs(C.momentum_x[id])); - min_l = std::min(min_l, std::abs(C.momentum_x[id])); + max_l = std::max(max_l, std::abs(C.momentum_x[id])); + min_l = std::min(min_l, std::abs(C.momentum_x[id])); } } } @@ -2570,29 +2567,28 @@ void Grid3D::Print_Grid_Stats(void) #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; + max_l = max_g; + min_l = min_g; #endif // MPI_CHOLLA - chprintf( - " abs(Momentum X) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", - mean_l, min_l, max_l); + chprintf(" abs(Momentum X) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", + mean_l, min_l, max_l); // y momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_y[id]); - max_l = std::max(max_l, std::abs(C.momentum_y[id])); - min_l = std::min(min_l, std::abs(C.momentum_y[id])); + max_l = std::max(max_l, std::abs(C.momentum_y[id])); + min_l = std::min(min_l, std::abs(C.momentum_y[id])); } } } @@ -2600,11 +2596,11 @@ void Grid3D::Print_Grid_Stats(void) #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; + max_l = max_g; + min_l = min_g; #endif // MPI_CHOLLA chprintf( " abs(Momentum Y) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", @@ -2613,16 +2609,16 @@ void Grid3D::Print_Grid_Stats(void) // z momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_z[id]); - max_l = std::max(max_l, std::abs(C.momentum_z[id])); - min_l = std::min(min_l, std::abs(C.momentum_z[id])); + max_l = std::max(max_l, std::abs(C.momentum_z[id])); + min_l = std::min(min_l, std::abs(C.momentum_z[id])); } } } @@ -2630,11 +2626,11 @@ void Grid3D::Print_Grid_Stats(void) #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; + max_l = max_g; + min_l = min_g; #endif // MPI_CHOLLA chprintf( " abs(Momentum Z) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", @@ -2642,16 +2638,16 @@ void Grid3D::Print_Grid_Stats(void) // Energy mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.Energy[id]; - max_l = std::max(max_l, C.Energy[id]); - min_l = std::min(min_l, C.Energy[id]); + max_l = std::max(max_l, C.Energy[id]); + min_l = std::min(min_l, C.Energy[id]); } } } @@ -2659,11 +2655,11 @@ void Grid3D::Print_Grid_Stats(void) #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; + max_l = max_g; + min_l = min_g; #endif // MPI_CHOLLA chprintf( " Energy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ]\n", @@ -2681,15 +2677,15 @@ void Grid3D::Print_Grid_Stats(void) for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); vx = C.momentum_x[id] / C.density[id]; vy = C.momentum_y[id] / C.density[id]; vz = C.momentum_z[id] / C.density[id]; gase = C.Energy[id] - 0.5 * C.density[id] * (vx * vx + vy * vy + vz * vz); mean_l += gase; - max_l = std::max(max_l, gase); - min_l = std::min(min_l, gase); + max_l = std::max(max_l, gase); + min_l = std::min(min_l, gase); temp = gase / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; temp_mean_l += temp; @@ -2732,11 +2728,11 @@ void Grid3D::Print_Grid_Stats(void) for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.GasEnergy[id]; - max_l = std::max(max_l, C.GasEnergy[id]); - min_l = std::min(min_l, C.GasEnergy[id]); + max_l = std::max(max_l, C.GasEnergy[id]); + min_l = std::min(min_l, C.GasEnergy[id]); temp = C.GasEnergy[id] / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; temp_mean_l += temp; @@ -2818,9 +2814,9 @@ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) #endif // DUST #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - #ifdef COSMOLOGY + #ifdef COSMOLOGY if (P.nfile > 0) { - #endif // COSMOLOGY + #endif // COSMOLOGY Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HI_density, "/HI_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HII_density, "/HII_density"); Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeI_density, "/HeI_density"); @@ -2830,27 +2826,26 @@ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) #ifdef GRACKLE_METALS Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.metal_density, "/metal_density"); #endif // GRACKLE_METALS - #ifdef COSMOLOGY + #ifdef COSMOLOGY } else { // Initialize the density field for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); - C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; - C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; - C.HeI_density[id] = INITIAL_FRACTION_HEI * C.density[id]; - C.HeII_density[id] = INITIAL_FRACTION_HEII * C.density[id]; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; + C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; + C.HeI_density[id] = INITIAL_FRACTION_HEI * C.density[id]; + C.HeII_density[id] = INITIAL_FRACTION_HEII * C.density[id]; C.HeIII_density[id] = INITIAL_FRACTION_HEIII * C.density[id]; - C.e_density[id] = INITIAL_FRACTION_ELECTRON * C.density[id]; + C.e_density[id] = INITIAL_FRACTION_ELECTRON * C.density[id]; } } } } - #endif // COSMOLOGY + #endif // COSMOLOGY #endif // COOLING_GRACKLE , CHEMISTRY_GPU - #endif // SCALAR // MHD only valid in 3D case diff --git a/src/io/io.h b/src/io/io.h index 8e3d5e66e..efee5c352 100644 --- a/src/io/io.h +++ b/src/io/io.h @@ -18,9 +18,8 @@ static inline bool Is_Root_Proc() #endif } -/*! \fn void Print_Stats(void) - * \brief Compute stats for a grid. */ -void Print_Stats(Grid3D &G); +/* Compute stats for a grid. */ +void Print_Stats(Grid3D& G); /* Write the data */ void Write_Data(Grid3D& G, struct Parameters P, int nfile); diff --git a/src/reconstruction/ppmp_cuda.cu b/src/reconstruction/ppmp_cuda.cu index 59efac774..6e0a13861 100644 --- a/src/reconstruction/ppmp_cuda.cu +++ b/src/reconstruction/ppmp_cuda.cu @@ -63,7 +63,7 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou // declare other variables Real del_q_imo, del_q_i, del_q_ipo; - //#ifndef VL + // #ifndef VL #ifdef CTU Real cs, cl, cr; // sound speed in cell i, and at left and right boundaries Real del_d, del_vx, del_vy, del_vz, del_p; // "slope" accross cell i From 48732b08fb8e5f2e8c58ae31d9b40a68f584829e Mon Sep 17 00:00:00 2001 From: Brant Robertson Date: Tue, 2 Jul 2024 13:56:35 -0700 Subject: [PATCH 61/62] More clang format tests --- src/io/io.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 236102fa8..0f4ccf8f7 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2497,9 +2497,9 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H } #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) -/*! \fn void Print_Stats(Grid3D& G) +/*! \fn void Print_Stats(Grid3D &G) * \brief Compute stats for a grid. */ -void Print_Stats(Grid3D& G) +void Print_Stats(Grid3D &G) { // Synchronize cudaMemcpy(G.C.density, G.C.device, G.H.n_fields * G.H.n_cells * sizeof(Real), cudaMemcpyDeviceToHost); @@ -2519,13 +2519,13 @@ void Grid3D::Print_Grid_Stats(void) // Density stats mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; // Do density first for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.density[id]; max_l = std::max(max_l, C.density[id]); From 04182aee08a6f4c8be967b3afd3655016d0328a8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 3 Jul 2024 21:13:51 +0000 Subject: [PATCH 62/62] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/analysis/io_analysis.cpp | 4 +- src/cosmology/cosmology.cpp | 4 +- src/cosmology/cosmology_functions.cpp | 2 - src/global/global.cpp | 5 +- src/io/io.cpp | 274 ++++++++++++-------------- src/reconstruction/ppmp_cuda.cu | 4 +- 6 files changed, 138 insertions(+), 155 deletions(-) diff --git a/src/analysis/io_analysis.cpp b/src/analysis/io_analysis.cpp index 88a94ccce..00efe714a 100644 --- a/src/analysis/io_analysis.cpp +++ b/src/analysis/io_analysis.cpp @@ -9,8 +9,8 @@ // #define OUTPUT_SKEWERS_TRANSMITTED_FLUX - using std::string; - using std::ifstream; +using std::ifstream; +using std::string; #ifdef OUTPUT_SKEWERS void Grid3D::Output_Skewers_File(struct Parameters *P) diff --git a/src/cosmology/cosmology.cpp b/src/cosmology/cosmology.cpp index a124303e5..2e63645c1 100644 --- a/src/cosmology/cosmology.cpp +++ b/src/cosmology/cosmology.cpp @@ -58,9 +58,7 @@ void Cosmology::Initialize(struct Parameters *P, Grav3D &Grav, Particles3D &Part Real dt_physical; // Advance a_t_sec until it matches current_a, and integrate time - while (a_t_sec < current_a) - { - + while (a_t_sec < current_a) { // Limit the scale a_t_sec factor to current_a if (a_t_sec + da_t_sec > current_a) { da_t_sec = current_a - a_t_sec; diff --git a/src/cosmology/cosmology_functions.cpp b/src/cosmology/cosmology_functions.cpp index f735258ce..76e7248a9 100644 --- a/src/cosmology/cosmology_functions.cpp +++ b/src/cosmology/cosmology_functions.cpp @@ -33,11 +33,9 @@ Real Cosmology::dtda_cosmo(Real da, Real a) return da / a_dot; } - /* Compute dt/da * da. dt/da is computed with a Runge-Kutta integration step */ Real Cosmology::Get_dt_from_da_rk(Real da, Real a) { - Real a3 = 0.3; Real a4 = 0.6; Real a5 = 1.0; diff --git a/src/global/global.cpp b/src/global/global.cpp index f5c758361..51df99757 100644 --- a/src/global/global.cpp +++ b/src/global/global.cpp @@ -103,7 +103,10 @@ char *Trim(char *s) // NOLINTNEXTLINE(cert-err58-cpp) // NOLINTNEXTLINE(*) -const std::set optionalParams = { "flag_delta","ddelta_dt","n_delta","Lz","Lx","phi","theta","delta","nzr","nxr","H0","Omega_M","Omega_L","Omega_R","Omega_K","w0","wa","Init_redshift","End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z" }; // NOLINT +const std::set optionalParams = { + "flag_delta", "ddelta_dt", "n_delta", "Lz", "Lx", "phi", "theta", "delta", + "nzr", "nxr", "H0", "Omega_M", "Omega_L", "Omega_R", "Omega_K", "w0", + "wa", "Init_redshift", "End_redshift", "tile_length", "n_proc_x", "n_proc_y", "n_proc_z"}; // NOLINT bool Old_Style_Parse_Param(const char *name, const char *value, struct Parameters *parms); diff --git a/src/io/io.cpp b/src/io/io.cpp index 0f4ccf8f7..7981dd54f 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -2496,7 +2496,7 @@ void Read_Grid_HDF5_Field_Magnetic(hid_t file_id, Real *dataset_buffer, Header H dataset_buffer, grid_buffer); } -#if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) + #if defined(PRINT_INITIAL_STATS) && defined(COSMOLOGY) /*! \fn void Print_Stats(Grid3D &G) * \brief Compute stats for a grid. */ void Print_Stats(Grid3D &G) @@ -2525,7 +2525,7 @@ void Grid3D::Print_Grid_Stats(void) for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.density[id]; max_l = std::max(max_l, C.density[id]); @@ -2535,27 +2535,26 @@ void Grid3D::Print_Grid_Stats(void) } mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); - #if MPI_CHOLLA + #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; - #endif // MPI_CHOLLA - chprintf("Density Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3] \n", - mean_l, min_l, max_l); + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf("Density Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3] \n", mean_l, min_l, max_l); // Momentum stats // x momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_x[id]); max_l = std::max(max_l, std::abs(C.momentum_x[id])); @@ -2565,26 +2564,24 @@ void Grid3D::Print_Grid_Stats(void) } mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); - #if MPI_CHOLLA + #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; - #endif // MPI_CHOLLA - chprintf(" abs(Momentum X) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", - mean_l, min_l, max_l); - + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf(" abs(Momentum X) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", mean_l, min_l, max_l); // y momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_y[id]); max_l = std::max(max_l, std::abs(C.momentum_y[id])); @@ -2594,27 +2591,24 @@ void Grid3D::Print_Grid_Stats(void) } mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); - #if MPI_CHOLLA + #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; - #endif // MPI_CHOLLA - chprintf( - " abs(Momentum Y) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", - mean_l, min_l, max_l); - + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf(" abs(Momentum Y) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", mean_l, min_l, max_l); // z momenta mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += std::abs(C.momentum_z[id]); max_l = std::max(max_l, std::abs(C.momentum_z[id])); @@ -2624,26 +2618,24 @@ void Grid3D::Print_Grid_Stats(void) } mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); - #if MPI_CHOLLA + #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; - #endif // MPI_CHOLLA - chprintf( - " abs(Momentum Z) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", - mean_l, min_l, max_l); + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf(" abs(Momentum Z) Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km s^-1] \n", mean_l, min_l, max_l); // Energy mean_l = 0; - min_l = 1e65; - max_l = -1; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.Energy[id]; max_l = std::max(max_l, C.Energy[id]); @@ -2653,120 +2645,112 @@ void Grid3D::Print_Grid_Stats(void) } mean_l /= ((H.nz_real) * (H.ny_real) * (H.nx_real)); - #if MPI_CHOLLA + #if MPI_CHOLLA mean_g = ReduceRealAvg(mean_l); - max_g = ReduceRealMax(max_l); - min_g = ReduceRealMin(min_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); mean_l = mean_g; - max_l = max_g; - min_l = min_g; - #endif // MPI_CHOLLA - chprintf( - " Energy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ]\n", - mean_l, min_l, max_l); + max_l = max_g; + min_l = min_g; + #endif // MPI_CHOLLA + chprintf(" Energy Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ]\n", mean_l, min_l, max_l); Real temp, temp_max_l, temp_min_l, temp_mean_l; Real temp_min_g, temp_max_g, temp_mean_g; Real gase, vx, vy, vz; temp_mean_l = 0; - temp_min_l = 1e65; - temp_max_l = -1; - mean_l = 0; - min_l = 1e65; - max_l = -1; + temp_min_l = 1e65; + temp_max_l = -1; + mean_l = 0; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); - vx = C.momentum_x[id] / C.density[id]; - vy = C.momentum_y[id] / C.density[id]; - vz = C.momentum_z[id] / C.density[id]; - gase = C.Energy[id] - 0.5 * C.density[id] * (vx * vx + vy * vy + vz * vz); + vx = C.momentum_x[id] / C.density[id]; + vy = C.momentum_y[id] / C.density[id]; + vz = C.momentum_z[id] / C.density[id]; + gase = C.Energy[id] - 0.5 * C.density[id] * (vx * vx + vy * vy + vz * vz); mean_l += gase; max_l = std::max(max_l, gase); min_l = std::min(min_l, gase); - temp = gase / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; + temp = gase / C.density[id] * (gama - 1) * MP / KB * 1e10; temp_mean_l += temp; temp_max_l = std::max(temp_max_l, temp); temp_min_l = std::min(temp_min_l, temp); } } } - mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); - temp_mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); - + mean_l /= (H.nz_real * H.ny_real * H.nx_real); + temp_mean_l /= (H.nz_real * H.ny_real * H.nx_real); - #if MPI_CHOLLA - mean_g = ReduceRealAvg( mean_l ); - max_g = ReduceRealMax( max_l ); - min_g = ReduceRealMin( min_l ); - mean_l = mean_g; - max_l = max_g; - min_l = min_g; - temp_mean_g = ReduceRealAvg( temp_mean_l ); - temp_max_g = ReduceRealMax( temp_max_l ); - temp_min_g = ReduceRealMin( temp_min_l ); + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + temp_mean_g = ReduceRealAvg(temp_mean_l); + temp_max_g = ReduceRealMax(temp_max_l); + temp_min_g = ReduceRealMin(temp_min_l); temp_mean_l = temp_mean_g; - temp_max_l = temp_max_g; - temp_min_l = temp_min_g; - #endif //MPI_CHOLLA + temp_max_l = temp_max_g; + temp_min_l = temp_min_g; + #endif // MPI_CHOLLA - chprintf( " GasEnergyCalc Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", - mean_l, min_l, max_l ); - chprintf( " TemperatureCalc Mean: %f Min: %f Max: %f [ K ] \n", - temp_mean_l, temp_min_l, temp_max_l ); + chprintf(" GasEnergyCalc Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", mean_l, min_l, max_l); + chprintf(" TemperatureCalc Mean: %f Min: %f Max: %f [ K ] \n", temp_mean_l, temp_min_l, temp_max_l); - #ifdef DE + #ifdef DE temp_mean_l = 0; - temp_min_l = 1e65; - temp_max_l = -1; - mean_l = 0; - min_l = 1e65; - max_l = -1; + temp_min_l = 1e65; + temp_max_l = -1; + mean_l = 0; + min_l = 1e65; + max_l = -1; for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); mean_l += C.GasEnergy[id]; max_l = std::max(max_l, C.GasEnergy[id]); min_l = std::min(min_l, C.GasEnergy[id]); - temp = C.GasEnergy[id] / C.density[id] * ( gama - 1 ) * MP / KB * 1e10 ; + temp = C.GasEnergy[id] / C.density[id] * (gama - 1) * MP / KB * 1e10; temp_mean_l += temp; temp_max_l = std::max(temp_max_l, temp); temp_min_l = std::min(temp_min_l, temp); } } } - mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); - temp_mean_l /= ( H.nz_real * H.ny_real * H.nx_real ); + mean_l /= (H.nz_real * H.ny_real * H.nx_real); + temp_mean_l /= (H.nz_real * H.ny_real * H.nx_real); - - #if MPI_CHOLLA - mean_g = ReduceRealAvg( mean_l ); - max_g = ReduceRealMax( max_l ); - min_g = ReduceRealMin( min_l ); - mean_l = mean_g; - max_l = max_g; - min_l = min_g; - temp_mean_g = ReduceRealAvg( temp_mean_l ); - temp_max_g = ReduceRealMax( temp_max_l ); - temp_min_g = ReduceRealMin( temp_min_l ); + #if MPI_CHOLLA + mean_g = ReduceRealAvg(mean_l); + max_g = ReduceRealMax(max_l); + min_g = ReduceRealMin(min_l); + mean_l = mean_g; + max_l = max_g; + min_l = min_g; + temp_mean_g = ReduceRealAvg(temp_mean_l); + temp_max_g = ReduceRealMax(temp_max_l); + temp_min_g = ReduceRealMin(temp_min_l); temp_mean_l = temp_mean_g; - temp_max_l = temp_max_g; - temp_min_l = temp_min_g; - #endif //MPI_CHOLLA - - chprintf( " GasEnergyDE Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", - mean_l, min_l, max_l ); - chprintf( " TemperatureDE Mean: %f Min: %f Max: %f [ K ] \n", - temp_mean_l, temp_min_l, temp_max_l ); - #endif // DE + temp_max_l = temp_max_g; + temp_min_l = temp_min_g; + #endif // MPI_CHOLLA + + chprintf(" GasEnergyDE Mean: %f Min: %f Max: %f [ h^2 Msun kpc^-3 km^2 s^-2 ] \n", mean_l, min_l, max_l); + chprintf(" TemperatureDE Mean: %f Min: %f Max: %f [ K ] \n", temp_mean_l, temp_min_l, temp_max_l); + #endif // DE } -#endif // PRINT_INITIAL_STATS and COSMOLOGY + #endif // PRINT_INITIAL_STATS and COSMOLOGY /*! \fn void Read_Grid_HDF5(hid_t file_id) * \brief Read in grid data from an hdf5 file. */ @@ -2813,40 +2797,40 @@ void Grid3D::Read_Grid_HDF5(hid_t file_id, struct Parameters P) Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.dust_density, "/dust_density"); #endif // DUST - #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) - #ifdef COSMOLOGY + #if defined(COOLING_GRACKLE) || defined(CHEMISTRY_GPU) + #ifdef COSMOLOGY if (P.nfile > 0) { - #endif // COSMOLOGY - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HI_density, "/HI_density"); - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HII_density, "/HII_density"); - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeI_density, "/HeI_density"); - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeII_density, "/HeII_density"); - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeIII_density, "/HeIII_density"); - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.e_density, "/e_density"); - #ifdef GRACKLE_METALS - Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.metal_density, "/metal_density"); - #endif // GRACKLE_METALS - #ifdef COSMOLOGY + #endif // COSMOLOGY + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HI_density, "/HI_density"); + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HII_density, "/HII_density"); + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeI_density, "/HeI_density"); + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeII_density, "/HeII_density"); + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.HeIII_density, "/HeIII_density"); + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.e_density, "/e_density"); + #ifdef GRACKLE_METALS + Read_Grid_HDF5_Field(file_id, dataset_buffer, H, C.metal_density, "/metal_density"); + #endif // GRACKLE_METALS + #ifdef COSMOLOGY } else { - // Initialize the density field + // Initialize the density field for (k = 0; k < H.nz_real; k++) { for (j = 0; j < H.ny_real; j++) { for (i = 0; i < H.nx_real; i++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; - buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); - C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; - C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; - C.HeI_density[id] = INITIAL_FRACTION_HEI * C.density[id]; - C.HeII_density[id] = INITIAL_FRACTION_HEII * C.density[id]; + id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + buf_id = k + j * (H.nz_real) + i * (H.nz_real) * (H.ny_real); + C.HI_density[id] = INITIAL_FRACTION_HI * C.density[id]; + C.HII_density[id] = INITIAL_FRACTION_HII * C.density[id]; + C.HeI_density[id] = INITIAL_FRACTION_HEI * C.density[id]; + C.HeII_density[id] = INITIAL_FRACTION_HEII * C.density[id]; C.HeIII_density[id] = INITIAL_FRACTION_HEIII * C.density[id]; - C.e_density[id] = INITIAL_FRACTION_ELECTRON * C.density[id]; + C.e_density[id] = INITIAL_FRACTION_ELECTRON * C.density[id]; } } } } - #endif // COSMOLOGY - #endif // COOLING_GRACKLE , CHEMISTRY_GPU - #endif // SCALAR + #endif // COSMOLOGY + #endif // COOLING_GRACKLE , CHEMISTRY_GPU + #endif // SCALAR // MHD only valid in 3D case if (H.nx > 1 && H.ny > 1 && H.nz > 1) { diff --git a/src/reconstruction/ppmp_cuda.cu b/src/reconstruction/ppmp_cuda.cu index 6e0a13861..e44a9ef42 100644 --- a/src/reconstruction/ppmp_cuda.cu +++ b/src/reconstruction/ppmp_cuda.cu @@ -83,7 +83,7 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou #ifdef DE Real ge_i, ge_imo, ge_ipo, ge_imt, ge_ipt, ge_L, ge_R, E_kin, E, dge; - // #ifndef VL + // #ifndef VL #ifdef CTU Real del_ge, ge_6, geL_0, geR_0; #endif // CTU @@ -92,7 +92,7 @@ __global__ void PPMP_cuda(Real *dev_conserved, Real *dev_bounds_L, Real *dev_bou #ifdef SCALAR Real scalar_i[NSCALARS], scalar_imo[NSCALARS], scalar_ipo[NSCALARS], scalar_imt[NSCALARS], scalar_ipt[NSCALARS]; Real scalar_L[NSCALARS], scalar_R[NSCALARS]; - // #ifndef VL + // #ifndef VL #ifdef CTU Real del_scalar[NSCALARS], scalar_6[NSCALARS], scalarL_0[NSCALARS], scalarR_0[NSCALARS]; #endif // CTU