Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add MHD Support to Projection and Rotated Projection Outputs #366

Merged
merged 9 commits into from
Jan 24, 2024
Merged

Add MHD Support to Projection and Rotated Projection Outputs #366

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
5fd419e
Add a function for computing the temperature from conserved variables
bcaddy Jan 22, 2024
ddd046e
Add MHD support to projection outputs
bcaddy Jan 23, 2024
80cd55e
Fix a DE bug in Dust_Kernel
bcaddy Jan 23, 2024
065cd79
Add MHD support to rotated projection outputs
bcaddy Jan 23, 2024
4226ddb
Fix test for new version of Calc_Temp_DE
bcaddy Jan 23, 2024
c84e4a2
Fix typo in Dust_Kernel, clarify comment
bcaddy Jan 23, 2024
f9bddfb
Replace `total_gas_energy` name with `gas_energy`
bcaddy Jan 24, 2024
ec6e34d
Refactor MHD support in `Dust_Kernel`
bcaddy Jan 24, 2024
495d453
Refactor temperature calculations in io.cpp
bcaddy Jan 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 22 additions & 35 deletions src/dust/dust_cuda.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -71,37 +66,33 @@ __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
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];
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);
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);
#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

// 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 tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) /
TIME_UNIT; // sputtering timescale, kyr (sim units)

Expand All @@ -121,10 +112,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
}
}

Expand Down
157 changes: 94 additions & 63 deletions src/io/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand All @@ -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;
Real const d = C.density[id];
dxy += d * 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 = 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
Txy += T * C.density[id] * H.dz;
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);
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;
}
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
Expand All @@ -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;
Real const d = C.density[id];
dxz += d * 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 = d * 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.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);
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;
}
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
Expand Down Expand Up @@ -1676,7 +1700,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;
Expand All @@ -1686,14 +1709,13 @@ 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, 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;
Real mu = 0.6;

srand(137); // initialize a random number
Expand Down Expand Up @@ -1740,11 +1762,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);
Expand All @@ -1769,33 +1794,39 @@ 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);
Real const 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);
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;
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;
}
}
}
Expand Down
38 changes: 35 additions & 3 deletions src/utils/hydro_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,43 @@ 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)
/*!
* \brief Compute the temperature when DE is turned on
*
* \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 gas_energy, Real const gamma, Real const n)
{
Real T = d * ge * (gamma - 1.0) * PRESSURE_UNIT / (n * KB);
return T;
return gas_energy * (gamma - 1.0) * PRESSURE_UNIT / (n * KB);
}
#endif // DE

Expand Down
Loading
Loading