From 5fd419e841a51f04a435a4a0c87235c06bb14f4a Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Mon, 22 Jan 2024 13:27:02 -0500 Subject: [PATCH 1/9] Add a function for computing the temperature from conserved variables --- src/utils/hydro_utilities.h | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index 24caff9f7..984a1f8ec 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -63,6 +63,30 @@ inline __host__ __device__ Real Calc_Temp(Real const &P, Real const &n) return T; } +/*! + * \brief Compute the temperature from the conserved variables + * + * \param[in] E The energy + * \param[in] d The density + * \param[in] mx The momentum in the X-direction + * \param[in] my The momentum in the Y-direction + * \param[in] mz The momentum in the Z-direction + * \param[in] gamma The adiabatic index + * \param[in] n The number density + * \param[in] magnetic_x The cell centered magnetic field in the X-direction + * \param[in] magnetic_y The cell centered magnetic field in the Y-direction + * \param[in] magnetic_z The cell centered magnetic field in the Z-direction + * \return Real The temperature of the gas in a cell + */ +inline __host__ __device__ Real Calc_Temp_Conserved(Real const E, Real const d, Real const mx, Real const my, + Real const mz, Real const gamma, Real const n, + Real const magnetic_x = 0.0, Real const magnetic_y = 0.0, + Real const magnetic_z = 0.0) +{ + Real const P = Calc_Pressure_Conserved(E, d, mx, my, mz, gamma, magnetic_x, magnetic_y, magnetic_z); + return Calc_Temp(P, n); +} + #ifdef DE inline __host__ __device__ Real Calc_Temp_DE(Real const &d, Real const &ge, Real const &gamma, Real const &n) { From ddd046e2398d9b38080f74f652b48864f4792de6 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 23 Jan 2024 13:56:16 -0500 Subject: [PATCH 2/9] Add MHD support to projection outputs --- src/io/io.cpp | 96 ++++++++++++++++++++++++++++++++------------------- 1 file changed, 60 insertions(+), 36 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index 9959267de..a246cb6b2 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -16,6 +16,8 @@ #include "../grid/grid3D.h" #include "../io/io.h" #include "../utils/cuda_utilities.h" +#include "../utils/hydro_utilities.h" +#include "../utils/mhd_utilities.h" #include "../utils/timing_functions.h" // provides ScopedTimer #ifdef MPI_CHOLLA #include "../mpi/mpi_routines.h" @@ -1526,18 +1528,16 @@ void Grid3D::Write_Grid_HDF5(hid_t file_id) * current simulation time. */ void Grid3D::Write_Projection_HDF5(hid_t file_id) { - int i, j, k, id, buf_id; hid_t dataset_id, dataspace_xy_id, dataspace_xz_id; Real *dataset_buffer_dxy, *dataset_buffer_dxz; Real *dataset_buffer_Txy, *dataset_buffer_Txz; herr_t status; - Real dxy, dxz, Txy, Txz, n, T; + Real dxy, dxz, Txy, Txz; #ifdef DUST Real dust_xy, dust_xz; Real *dataset_buffer_dust_xy, *dataset_buffer_dust_xz; #endif - n = T = 0; Real mu = 0.6; // 3D @@ -1563,37 +1563,51 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) dataspace_xz_id = H5Screate_simple(2, dims, NULL); // Copy the xy density and temperature projections to the memory buffer - for (j = 0; j < H.ny_real; j++) { - for (i = 0; i < H.nx_real; i++) { + for (int j = 0; j < H.ny_real; j++) { + for (int i = 0; i < H.nx_real; i++) { dxy = 0; Txy = 0; #ifdef DUST dust_xy = 0; #endif // for each xy element, sum over the z column - for (k = 0; k < H.nz_real; k++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (int k = 0; k < H.nz_real; k++) { + int const xid = i + H.n_ghost; + int const yid = j + H.n_ghost; + int const zid = k + H.n_ghost; + int const id = cuda_utilities::compute1DIndex(xid, yid, zid, H.nx, H.ny); + // sum density dxy += C.density[id] * H.dz; #ifdef DUST dust_xy += C.dust_density[id] * H.dz; #endif // calculate number density - n = C.density[id] * DENSITY_UNIT / (mu * MP); + Real const n = C.density[id] * DENSITY_UNIT / (mu * MP); + // calculate temperature - #ifndef DE - Real mx = C.momentum_x[id]; - Real my = C.momentum_y[id]; - Real mz = C.momentum_z[id]; - Real E = C.Energy[id]; - T = (E - 0.5 * (mx * mx + my * my + mz * mz) / C.density[id]) * (gama - 1.0) * PRESSURE_UNIT / (n * KB); - #endif #ifdef DE - T = C.GasEnergy[id] * PRESSURE_UNIT * (gama - 1.0) / (n * KB); - #endif - Txy += T * C.density[id] * H.dz; + Real const T = hydro_utilities::Calc_Temp_DE(C.density[id], C.GasEnergy[id], gama, n); + #else // DE is not defined + Real const mx = C.momentum_x[id]; + Real const my = C.momentum_y[id]; + Real const mz = C.momentum_z[id]; + Real const E = C.Energy[id]; + + #ifdef MHD + auto const [magnetic_x, magnetic_y, magnetic_z] = + mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); + #else // MHD is not defined + Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; + #endif // MHD + + Real const T = hydro_utilities::Calc_Temp_Conserved(E, C.density[id], mx, my, mz, gama, n, magnetic_x, + magnetic_y, magnetic_z); + #endif // DE + + Txy += T * H.dz; } - buf_id = j + i * H.ny_real; + int const buf_id = j + i * H.ny_real; dataset_buffer_dxy[buf_id] = dxy; dataset_buffer_Txy[buf_id] = Txy; #ifdef DUST @@ -1603,37 +1617,47 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) } // Copy the xz density and temperature projections to the memory buffer - for (k = 0; k < H.nz_real; k++) { - for (i = 0; i < H.nx_real; i++) { + for (int k = 0; k < H.nz_real; k++) { + for (int i = 0; i < H.nx_real; i++) { dxz = 0; Txz = 0; #ifdef DUST dust_xz = 0; #endif // for each xz element, sum over the y column - for (j = 0; j < H.ny_real; j++) { - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + for (int j = 0; j < H.ny_real; j++) { + int const xid = i + H.n_ghost; + int const yid = j + H.n_ghost; + int const zid = k + H.n_ghost; + int const id = cuda_utilities::compute1DIndex(xid, yid, zid, H.nx, H.ny); // sum density dxz += C.density[id] * H.dy; #ifdef DUST dust_xz += C.dust_density[id] * H.dy; #endif // calculate number density - n = C.density[id] * DENSITY_UNIT / (mu * MP); - // calculate temperature - #ifndef DE - Real mx = C.momentum_x[id]; - Real my = C.momentum_y[id]; - Real mz = C.momentum_z[id]; - Real E = C.Energy[id]; - T = (E - 0.5 * (mx * mx + my * my + mz * mz) / C.density[id]) * (gama - 1.0) * PRESSURE_UNIT / (n * KB); - #endif + Real const n = C.density[id] * DENSITY_UNIT / (mu * MP); #ifdef DE - T = C.GasEnergy[id] * PRESSURE_UNIT * (gama - 1.0) / (n * KB); - #endif - Txz += T * C.density[id] * H.dy; + Real const T = hydro_utilities::Calc_Temp_DE(C.density[id], C.GasEnergy[id], gama, n); + #else // DE is not defined + Real const mx = C.momentum_x[id]; + Real const my = C.momentum_y[id]; + Real const mz = C.momentum_z[id]; + Real const E = C.Energy[id]; + + #ifdef MHD + auto const [magnetic_x, magnetic_y, magnetic_z] = + mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); + #else // MHD is not defined + Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; + #endif // MHD + + Real const T = hydro_utilities::Calc_Temp_Conserved(E, C.density[id], mx, my, mz, gama, n, magnetic_x, + magnetic_y, magnetic_z); + #endif // DE + Txz += T * H.dy; } - buf_id = k + i * H.nz_real; + int const buf_id = k + i * H.nz_real; dataset_buffer_dxz[buf_id] = dxz; dataset_buffer_Txz[buf_id] = Txz; #ifdef DUST From 80cd55e5b3d641b8d209be8b87db7c278705c508 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 23 Jan 2024 16:17:51 -0500 Subject: [PATCH 3/9] Fix a DE bug in Dust_Kernel The gas energy was being loaded then checked for correctness against a variable that didn't exist. The gas energy was also being updated when it shouldn't be. Fixed those bugs and added MHD support to the temperature calculations. Also, updated the hydro_utilities::Calc_Temp_DE function to use the total, not specific, gas energy and added documentation for clarity --- src/dust/dust_cuda.cu | 59 ++++++++++++++----------------------- src/io/io.cpp | 26 ++++++++-------- src/utils/hydro_utilities.h | 13 ++++++-- 3 files changed, 46 insertions(+), 52 deletions(-) diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu index 46273ef03..ced77858c 100644 --- a/src/dust/dust_cuda.cu +++ b/src/dust/dust_cuda.cu @@ -52,14 +52,9 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g int id_x = id - id_z * nx * ny - id_y * nx; // define physics variables - Real density_gas, density_dust; // fluid mass densities - Real number_density; // gas number density - Real mu = 0.6; // mean molecular weight - Real temperature, energy, pressure; // temperature, energy, pressure - Real velocity_x, velocity_y, velocity_z; // velocities - #ifdef DE - Real energy_gas; - #endif // DE + Real density_gas, density_dust; // fluid mass densities + Real number_density; // gas number density + Real mu = 0.6; // mean molecular weight // define integration variables Real dd_dt; // instantaneous rate of change in dust density @@ -71,36 +66,30 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g // get conserved quanitites density_gas = dev_conserved[id + n_cells * grid_enum::density]; density_dust = dev_conserved[id + n_cells * grid_enum::dust_density]; - energy = dev_conserved[id + n_cells * grid_enum::Energy]; // convert mass density to number density number_density = density_gas * DENSITY_UNIT / (mu * MP); - if (energy < 0.0 || energy != energy) { - return; - } - - // get conserved quanitites - velocity_x = dev_conserved[id + n_cells * grid_enum::momentum_x] / density_gas; - velocity_y = dev_conserved[id + n_cells * grid_enum::momentum_y] / density_gas; - velocity_z = dev_conserved[id + n_cells * grid_enum::momentum_z] / density_gas; + // Compute the temperature #ifdef DE - energy_gas = dev_conserved[id + n_cells * grid_enum::GasEnergy] / density_gas; - energy_gas = fmax(ge, (Real)TINY_NUMBER); - #endif // DE - - // calculate physical quantities - pressure = hydro_utilities::Calc_Pressure_Primitive(energy, density_gas, velocity_x, velocity_y, velocity_z, gamma); - - Real temperature_init; - temperature_init = hydro_utilities::Calc_Temp(pressure, number_density); - - #ifdef DE - temperature_init = hydro_utilities::Calc_Temp_DE(density_gas, energy_gas, gamma, number_density); - #endif // DE - - // if dual energy is turned on use temp from total internal energy - temperature = temperature_init; + Real const total_gas_energy = dev_conserved[id + n_cells * grid_enum::GasEnergy]; + Real const temperature = hydro_utilities::Calc_Temp_DE(total_gas_energy, gamma, number_density); + #else // DE is not enabled + Real const energy = dev_conserved[id + n_cells * grid_enum::Energy]; + Real const momentum_x = dev_conserved[id + n_cells * grid_enum::momentum_x]; + Real const momentum_y = dev_conserved[id + n_cells * grid_enum::momentum_y]; + Real const momentum_z = dev_conserved[id + n_cells * grid_enum::momentum_z]; + + #ifdef MHD + auto const [magnetic_x, magnetic_y, magnetic_z] = + mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); + #else // MHD is not defined + Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; + #endif // MHD + + Real const temperature = hydro_utilities::Calc_Temp_Conserved( + energy, density, momentum_x, momentum_y, momentum_z, gamma, number_density, magnetic_x, magnetic_y, magnetic_z); + #endif // DE Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) / TIME_UNIT; // sputtering timescale, kyr (sim units) @@ -121,10 +110,6 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g density_dust += dd; dev_conserved[id + n_cells * grid_enum::dust_density] = density_dust; - - #ifdef DE - dev_conserved[id + n_cells * grid_enum::GasEnergy] = density_dust * energy_gas; - #endif } } diff --git a/src/io/io.cpp b/src/io/io.cpp index a246cb6b2..a16267867 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -1578,16 +1578,17 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) int const id = cuda_utilities::compute1DIndex(xid, yid, zid, H.nx, H.ny); // sum density - dxy += C.density[id] * H.dz; + Real const d = C.density[id]; + dxy += d * H.dz; #ifdef DUST dust_xy += C.dust_density[id] * H.dz; #endif // calculate number density - Real const n = C.density[id] * DENSITY_UNIT / (mu * MP); + Real const n = d * DENSITY_UNIT / (mu * MP); // calculate temperature #ifdef DE - Real const T = hydro_utilities::Calc_Temp_DE(C.density[id], C.GasEnergy[id], gama, n); + Real const T = hydro_utilities::Calc_Temp_DE(C.GasEnergy[id], gama, n); #else // DE is not defined Real const mx = C.momentum_x[id]; Real const my = C.momentum_y[id]; @@ -1601,11 +1602,11 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; #endif // MHD - Real const T = hydro_utilities::Calc_Temp_Conserved(E, C.density[id], mx, my, mz, gama, n, magnetic_x, - magnetic_y, magnetic_z); + Real const T = + hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); #endif // DE - Txy += T * H.dz; + Txy += T * d * H.dz; } int const buf_id = j + i * H.ny_real; dataset_buffer_dxy[buf_id] = dxy; @@ -1631,14 +1632,15 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) int const zid = k + H.n_ghost; int const id = cuda_utilities::compute1DIndex(xid, yid, zid, H.nx, H.ny); // sum density - dxz += C.density[id] * H.dy; + Real const d = C.density[id]; + dxz += d * H.dy; #ifdef DUST dust_xz += C.dust_density[id] * H.dy; #endif // calculate number density - Real const n = C.density[id] * DENSITY_UNIT / (mu * MP); + Real const n = d * DENSITY_UNIT / (mu * MP); #ifdef DE - Real const T = hydro_utilities::Calc_Temp_DE(C.density[id], C.GasEnergy[id], gama, n); + Real const T = hydro_utilities::Calc_Temp_DE(C.GasEnergy[id], gama, n); #else // DE is not defined Real const mx = C.momentum_x[id]; Real const my = C.momentum_y[id]; @@ -1652,10 +1654,10 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; #endif // MHD - Real const T = hydro_utilities::Calc_Temp_Conserved(E, C.density[id], mx, my, mz, gama, n, magnetic_x, - magnetic_y, magnetic_z); + Real const T = + hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); #endif // DE - Txz += T * H.dy; + Txz += T * d * H.dy; } int const buf_id = k + i * H.nz_real; dataset_buffer_dxz[buf_id] = dxz; diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index 984a1f8ec..e49554366 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -88,10 +88,17 @@ inline __host__ __device__ Real Calc_Temp_Conserved(Real const E, Real const d, } #ifdef DE -inline __host__ __device__ Real Calc_Temp_DE(Real const &d, Real const &ge, Real const &gamma, Real const &n) +/*! + * \brief Compute the temperature when DE is turned on + * + * \param total_gas_energy The total gas energy in the cell. This is the value stored in the grid + * \param gamma The adiabatic index + * \param n The number density + * \return Real The temperature + */ +inline __host__ __device__ Real Calc_Temp_DE(Real const &total_gas_energy, Real const &gamma, Real const &n) { - Real T = d * ge * (gamma - 1.0) * PRESSURE_UNIT / (n * KB); - return T; + return total_gas_energy * (gamma - 1.0) * PRESSURE_UNIT / (n * KB); } #endif // DE From 065cd79b691ad967614298bbdcd59f238c5fa711 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 23 Jan 2024 16:35:10 -0500 Subject: [PATCH 4/9] Add MHD support to rotated projection outputs --- src/io/io.cpp | 57 +++++++++++++++++++++++++++++---------------------- 1 file changed, 33 insertions(+), 24 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index a16267867..f93a30360 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -1702,7 +1702,6 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) * time. */ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) { - int i, j, k, id, buf_id; hid_t dataset_id, dataspace_xzr_id; Real *dataset_buffer_dxzr; Real *dataset_buffer_Txzr; @@ -1712,14 +1711,14 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) herr_t status; Real dxy, dxz, Txy, Txz; - Real d, n, T, vx, vy, vz; + Real d, n, vx, vy, vz; Real x, y, z; // cell positions Real xp, yp, zp; // rotated positions Real alpha, beta; // projected positions int ix, iz; // projected index positions - n = T = 0; + n = 0; Real mu = 0.6; srand(137); // initialize a random number @@ -1766,11 +1765,14 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) dataspace_xzr_id = H5Screate_simple(2, dims, NULL); // Copy the xz rotated projection to the memory buffer - for (k = 0; k < H.nz_real; k++) { - for (i = 0; i < H.nx_real; i++) { - for (j = 0; j < H.ny_real; j++) { + for (int k = 0; k < H.nz_real; k++) { + for (int i = 0; i < H.nx_real; i++) { + for (int j = 0; j < H.ny_real; j++) { // get cell index - id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny; + int const xid = i + H.n_ghost; + int const yid = j + H.n_ghost; + int const zid = k + H.n_ghost; + int const id = cuda_utilities::compute1DIndex(xid, yid, zid, H.nx, H.ny); // get cell positions Get_Position(i + H.n_ghost, j + H.n_ghost, k + H.n_ghost, &x, &y, &z); @@ -1795,33 +1797,40 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) #endif if ((ix >= 0) && (ix < nx_dset) && (iz >= 0) && (iz < nz_dset)) { - buf_id = iz + ix * nz_dset; - d = C.density[id]; + int const buf_id = iz + ix * nz_dset; + d = C.density[id]; // project density dataset_buffer_dxzr[buf_id] += d * H.dy; // calculate number density n = d * DENSITY_UNIT / (mu * MP); + // calculate temperature - #ifndef DE - Real mx = C.momentum_x[id]; - Real my = C.momentum_y[id]; - Real mz = C.momentum_z[id]; - Real E = C.Energy[id]; - T = (E - 0.5 * (mx * mx + my * my + mz * mz) / C.density[id]) * (gama - 1.0) * PRESSURE_UNIT / (n * KB); - #endif #ifdef DE - T = C.GasEnergy[id] * PRESSURE_UNIT * (gama - 1.0) / (n * KB); - #endif + Real const T = hydro_utilities::Calc_Temp_DE(C.GasEnergy[id], gama, n); + #else // DE is not defined + Real const mx = C.momentum_x[id]; + Real const my = C.momentum_y[id]; + Real const mz = C.momentum_z[id]; + Real const E = C.Energy[id]; + + #ifdef MHD + auto const [magnetic_x, magnetic_y, magnetic_z] = + mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); + #else // MHD is not defined + Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; + #endif // MHD + + Real const T = + hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); + #endif // DE + Txz = T * d * H.dy; dataset_buffer_Txzr[buf_id] += Txz; // compute velocities - vx = C.momentum_x[id]; - dataset_buffer_vxxzr[buf_id] += vx * H.dy; - vy = C.momentum_y[id]; - dataset_buffer_vyxzr[buf_id] += vy * H.dy; - vz = C.momentum_z[id]; - dataset_buffer_vzxzr[buf_id] += vz * H.dy; + dataset_buffer_vxxzr[buf_id] += C.momentum_x[id] * H.dy; + dataset_buffer_vyxzr[buf_id] += C.momentum_y[id] * H.dy; + dataset_buffer_vzxzr[buf_id] += C.momentum_z[id] * H.dy; } } } From 4226ddba591673b6da7dab165b52007b45e5241a Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 23 Jan 2024 16:37:54 -0500 Subject: [PATCH 5/9] Fix test for new version of Calc_Temp_DE --- src/utils/hydro_utilities_tests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/hydro_utilities_tests.cpp b/src/utils/hydro_utilities_tests.cpp index b200ddd8c..fe0cbe9e6 100644 --- a/src/utils/hydro_utilities_tests.cpp +++ b/src/utils/hydro_utilities_tests.cpp @@ -145,7 +145,7 @@ TEST(tHYDROHydroUtilsCalcTempDE, CorrectInputExpectCorrectOutput) for (size_t i = 0; i < parameters.names.size(); i++) { Real test_Ts = - hydro_utilities::Calc_Temp_DE(parameters.d.at(i), parameters.ge.at(i), parameters.gamma, parameters.n.at(i)); + hydro_utilities::Calc_Temp_DE(parameters.d.at(i) * parameters.ge.at(i), parameters.gamma, parameters.n.at(i)); testing_utilities::Check_Results(fiducial_Ts.at(i), test_Ts, parameters.names.at(i)); } From c84e4a292b4ab49ae0161773dbe1da47558916c5 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Tue, 23 Jan 2024 16:46:32 -0500 Subject: [PATCH 6/9] Fix typo in Dust_Kernel, clarify comment --- src/dust/dust_cuda.cu | 5 +++-- src/utils/hydro_utilities.h | 7 ++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu index ced77858c..c38180ab5 100644 --- a/src/dust/dust_cuda.cu +++ b/src/dust/dust_cuda.cu @@ -87,8 +87,9 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; #endif // MHD - Real const temperature = hydro_utilities::Calc_Temp_Conserved( - energy, density, momentum_x, momentum_y, momentum_z, gamma, number_density, magnetic_x, magnetic_y, magnetic_z); + Real const temperature = + hydro_utilities::Calc_Temp_Conserved(energy, density_gas, momentum_x, momentum_y, momentum_z, gamma, + number_density, magnetic_x, magnetic_y, magnetic_z); #endif // DE Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) / diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index e49554366..fab5aece1 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -91,9 +91,10 @@ inline __host__ __device__ Real Calc_Temp_Conserved(Real const E, Real const d, /*! * \brief Compute the temperature when DE is turned on * - * \param total_gas_energy The total gas energy in the cell. This is the value stored in the grid - * \param gamma The adiabatic index - * \param n The number density + * \param[in] total_gas_energy The total gas energy in the cell. This is the value stored in the grid at + * grid_enum::GasEnergy + * \param[in] gamma The adiabatic index + * \param[in] n The number density * \return Real The temperature */ inline __host__ __device__ Real Calc_Temp_DE(Real const &total_gas_energy, Real const &gamma, Real const &n) From f9bddfb684b9c694fdffbf7a0523a50e7ad0729d Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Wed, 24 Jan 2024 11:04:42 -0500 Subject: [PATCH 7/9] Replace `total_gas_energy` name with `gas_energy` gas_energy is the total gas energy and specific_gas_energy is the specific version --- src/dust/dust_cuda.cu | 4 ++-- src/utils/hydro_utilities.h | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu index c38180ab5..4dc05648c 100644 --- a/src/dust/dust_cuda.cu +++ b/src/dust/dust_cuda.cu @@ -72,8 +72,8 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g // Compute the temperature #ifdef DE - Real const total_gas_energy = dev_conserved[id + n_cells * grid_enum::GasEnergy]; - Real const temperature = hydro_utilities::Calc_Temp_DE(total_gas_energy, gamma, number_density); + Real const gas_energy = dev_conserved[id + n_cells * grid_enum::GasEnergy]; + Real const temperature = hydro_utilities::Calc_Temp_DE(gas_energy, gamma, number_density); #else // DE is not enabled Real const energy = dev_conserved[id + n_cells * grid_enum::Energy]; Real const momentum_x = dev_conserved[id + n_cells * grid_enum::momentum_x]; diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h index fab5aece1..1a464e899 100644 --- a/src/utils/hydro_utilities.h +++ b/src/utils/hydro_utilities.h @@ -91,15 +91,15 @@ inline __host__ __device__ Real Calc_Temp_Conserved(Real const E, Real const d, /*! * \brief Compute the temperature when DE is turned on * - * \param[in] total_gas_energy The total gas energy in the cell. This is the value stored in the grid at + * \param[in] gas_energy The total gas energy in the cell. This is the value stored in the grid at * grid_enum::GasEnergy * \param[in] gamma The adiabatic index * \param[in] n The number density * \return Real The temperature */ -inline __host__ __device__ Real Calc_Temp_DE(Real const &total_gas_energy, Real const &gamma, Real const &n) +inline __host__ __device__ Real Calc_Temp_DE(Real const gas_energy, Real const gamma, Real const n) { - return total_gas_energy * (gamma - 1.0) * PRESSURE_UNIT / (n * KB); + return gas_energy * (gamma - 1.0) * PRESSURE_UNIT / (n * KB); } #endif // DE From ec6e34d378e1683e13b000ccd3121be0462cb1a1 Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Wed, 24 Jan 2024 11:06:27 -0500 Subject: [PATCH 8/9] Refactor MHD support in `Dust_Kernel` Now you only use magnetic variables when MHD is enabled --- src/dust/dust_cuda.cu | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu index 4dc05648c..69cc83eec 100644 --- a/src/dust/dust_cuda.cu +++ b/src/dust/dust_cuda.cu @@ -83,14 +83,15 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g #ifdef MHD auto const [magnetic_x, magnetic_y, magnetic_z] = mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); - #else // MHD is not defined - Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; - #endif // MHD - Real const temperature = hydro_utilities::Calc_Temp_Conserved(energy, density_gas, momentum_x, momentum_y, momentum_z, gamma, number_density, magnetic_x, magnetic_y, magnetic_z); - #endif // DE + #else // MHD is not defined + Real const temperature = hydro_utilities::Calc_Temp_Conserved(energy, density_gas, momentum_x, momentum_y, + momentum_z, gamma, number_density); + #endif // MHD + + #endif // DE Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) / TIME_UNIT; // sputtering timescale, kyr (sim units) From 495d453e0e2e3fbad687988b2b1d7d35d10d867a Mon Sep 17 00:00:00 2001 From: Bob Caddy Date: Wed, 24 Jan 2024 12:53:30 -0500 Subject: [PATCH 9/9] Refactor temperature calculations in io.cpp --- src/io/io.cpp | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/io/io.cpp b/src/io/io.cpp index f93a30360..51303c72d 100644 --- a/src/io/io.cpp +++ b/src/io/io.cpp @@ -1598,12 +1598,11 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) #ifdef MHD auto const [magnetic_x, magnetic_y, magnetic_z] = mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); - #else // MHD is not defined - Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; - #endif // MHD - Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); + #else // MHD is not defined + Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n); + #endif // MHD #endif // DE Txy += T * d * H.dz; @@ -1650,12 +1649,11 @@ void Grid3D::Write_Projection_HDF5(hid_t file_id) #ifdef MHD auto const [magnetic_x, magnetic_y, magnetic_z] = mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); - #else // MHD is not defined - Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; - #endif // MHD - Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); + #else // MHD is not defined + Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n); + #endif // MHD #endif // DE Txz += T * d * H.dy; } @@ -1711,14 +1709,13 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) herr_t status; Real dxy, dxz, Txy, Txz; - Real d, n, vx, vy, vz; + Real d, vx, vy, vz; Real x, y, z; // cell positions Real xp, yp, zp; // rotated positions Real alpha, beta; // projected positions int ix, iz; // projected index positions - n = 0; Real mu = 0.6; srand(137); // initialize a random number @@ -1802,7 +1799,7 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) // project density dataset_buffer_dxzr[buf_id] += d * H.dy; // calculate number density - n = d * DENSITY_UNIT / (mu * MP); + Real const n = d * DENSITY_UNIT / (mu * MP); // calculate temperature #ifdef DE @@ -1816,12 +1813,11 @@ void Grid3D::Write_Rotated_Projection_HDF5(hid_t file_id) #ifdef MHD auto const [magnetic_x, magnetic_y, magnetic_z] = mhd::utils::cellCenteredMagneticFields(C.host, id, xid, yid, zid, H.n_cells, H.nx, H.ny); - #else // MHD is not defined - Real const magnetic_x = 0.0, magnetic_y = 0.0, magnetic_z = 0.0; - #endif // MHD - Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n, magnetic_x, magnetic_y, magnetic_z); + #else // MHD is not defined + Real const T = hydro_utilities::Calc_Temp_Conserved(E, d, mx, my, mz, gama, n); + #endif // MHD #endif // DE Txz = T * d * H.dy;