Skip to content

Commit

Permalink
Merge pull request #366 from bcaddy/dev-mhdProjections
Browse files Browse the repository at this point in the history 
Add MHD Support to Projection and Rotated Projection Outputs
  • Loading branch information
@evaneschneider
evaneschneider authored Jan 24, 2024
2 parents db10522 + 495d453 commit 70341bd
Show file tree
Hide file tree
Showing 4 changed files with 152 additions and 102 deletions.
57 changes: 22 additions & 35 deletions src/dust/dust_cuda.cu
View file
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
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
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

0 comments on commit 70341bd

Please sign in to comment.