mhd_utilities.h

Go to the documentation of this file.

 
#pragma once
 
// STL Includes
#include <vector>
 
// External Includes
 
// Local Includes
#include "../global/global.h"
#include "../global/global_cuda.h"
#include "../grid/grid3D.h"
#include "../utils/basic_structs.h"
#include "../utils/cuda_utilities.h"
#include "../utils/gpu.hpp"
#include "../utils/math_utilities.h"
 
namespace mhd::utils
{
namespace internal
{
// =====================================================================
inline __host__ __device__ Real _magnetosonicSpeed(Real const &density, Real const &gasPressure, Real const &magneticX,
                                                   Real const &magneticY, Real const &magneticZ, Real const &gamma,
                                                   Real const &waveChoice)
{
  // Compute the sound speed
  Real bXSquared = magneticX * magneticX;
  Real bSquared  = bXSquared + ((magneticY * magneticY) + (magneticZ * magneticZ));
 
  Real term1 = gamma * gasPressure + bSquared;
 
  Real term2 = (term1 * term1) - 4. * gamma * gasPressure * bXSquared;
  term2      = sqrt(term2);
 
  return sqrt((term1 + waveChoice * term2) / (2.0 * fmax(density, TINY_NUMBER)));
}
// =====================================================================
}  // namespace internal
 
// =========================================================================
inline __host__ __device__ Real computeMagneticEnergy(Real const &magneticX, Real const &magneticY,
                                                      Real const &magneticZ)
{
  return 0.5 * math_utils::SquareMagnitude(magneticX, magneticY, magneticZ);
}
// =========================================================================
 
// =========================================================================
inline __host__ __device__ Real computeThermalEnergy(Real const &energyTot, Real const &density, Real const &momentumX,
                                                     Real const &momentumY, Real const &momentumZ,
                                                     Real const &magneticX, Real const &magneticY,
                                                     Real const &magneticZ, Real const &gamma)
{
  return energyTot - 0.5 * math_utils::SquareMagnitude(momentumX, momentumY, momentumZ) / fmax(density, TINY_NUMBER) -
         computeMagneticEnergy(magneticX, magneticY, magneticZ);
}
// =========================================================================
 
// =========================================================================
inline __host__ __device__ Real computeTotalPressure(Real const &gasPressure, Real const &magneticX,
                                                     Real const &magneticY, Real const &magneticZ)
{
  Real pTot = gasPressure + computeMagneticEnergy(magneticX, magneticY, magneticZ);
 
  return fmax(pTot, TINY_NUMBER);
}
inline __host__ __device__ Real computeTotalPressure(Real const &gasPressure,
                                                     hydro_utilities::VectorXYZ<Real> const &magnetic)
{
  return computeTotalPressure(gasPressure, magnetic.x(), magnetic.y(), magnetic.z());
}
// =========================================================================
 
// =========================================================================
inline __host__ __device__ Real fastMagnetosonicSpeed(Real const &density, Real const &pressure, Real const &magneticX,
                                                      Real const &magneticY, Real const &magneticZ, Real const &gamma)
{
  // Compute the sound speed
  return mhd::utils::internal::_magnetosonicSpeed(density, pressure, magneticX, magneticY, magneticZ, gamma, 1.0);
}
// =========================================================================
 
// =========================================================================
inline __host__ __device__ Real slowMagnetosonicSpeed(Real const &density, Real const &pressure, Real const &magneticX,
                                                      Real const &magneticY, Real const &magneticZ, Real const &gamma)
{
  // Compute the sound speed
  return mhd::utils::internal::_magnetosonicSpeed(density, pressure, magneticX, magneticY, magneticZ, gamma, -1.0);
}
// =========================================================================
 
// =========================================================================
inline __host__ __device__ Real alfvenSpeed(Real const &magneticX, Real const &density)
{
  // Compute the Alfven wave speed
  return fabs(magneticX) / sqrt(fmax(density, TINY_NUMBER));
}
// =========================================================================
 
// =========================================================================
#ifdef MHD
inline __host__ __device__ hydro_utilities::VectorXYZ<Real> cellCenteredMagneticFields(
    Real const *dev_conserved, size_t const &id, size_t const &xid, size_t const &yid, size_t const &zid,
    size_t const &n_cells, size_t const &nx, size_t const &ny)
{
  // Ternary operator to check that no values outside of the magnetic field
  // arrays are loaded. If the cell is on the edge that doesn't have magnetic
  // fields on both sides then instead set the centered magnetic field to be
  // equal to the magnetic field of the closest edge.
  Real avgBx = (xid > 0) ?
                         /*if true*/ 0.5 * (dev_conserved[(grid_enum::magnetic_x)*n_cells + id] +
                                            dev_conserved[(grid_enum::magnetic_x)*n_cells +
                                                          cuda_utilities::compute1DIndex(xid - 1, yid, zid, nx, ny)])
                         :
                         /*if false*/ dev_conserved[(grid_enum::magnetic_x)*n_cells + id];
  Real avgBy = (yid > 0) ?
                         /*if true*/ 0.5 * (dev_conserved[(grid_enum::magnetic_y)*n_cells + id] +
                                            dev_conserved[(grid_enum::magnetic_y)*n_cells +
                                                          cuda_utilities::compute1DIndex(xid, yid - 1, zid, nx, ny)])
                         :
                         /*if false*/ dev_conserved[(grid_enum::magnetic_y)*n_cells + id];
  Real avgBz = (zid > 0) ?
                         /*if true*/ 0.5 * (dev_conserved[(grid_enum::magnetic_z)*n_cells + id] +
                                            dev_conserved[(grid_enum::magnetic_z)*n_cells +
                                                          cuda_utilities::compute1DIndex(xid, yid, zid - 1, nx, ny)])
                         :
                         /*if false*/ dev_conserved[(grid_enum::magnetic_z)*n_cells + id];
 
  return {avgBx, avgBy, avgBz};
}
// =========================================================================
 
// =========================================================================
void Init_Magnetic_Field_With_Vector_Potential(Header const &H, Grid3D::Conserved const &C,
                                               std::vector<Real> const &vectorPotential);
// =========================================================================
#endif  // MHD
}  // end namespace mhd::utils