dual_energy.h

Go to the documentation of this file.

 
#ifndef DUAL_ENERGY_H
#define DUAL_ENERGY_H
 
#include "../global/global.h"
#include "../grid/grid_enum.h"
#include "../hydro/hydro_cuda.h"
#include "../utils/basic_structs.h"
 
namespace dual_energy
{
 
template <int NDim>
__global__ void Select_Internal_Energy(Real *dev_conserved, hydro_utilities::VectorXYZ<int> grid_shape, int n_ghost,
                                       int n_fields)
{
  static_assert((NDim == 1) || (NDim == 2) || (NDim == 3), "NDim must be 1, 2, or 3");
 
#ifndef DE
  printf("WARNING: this function isn't usable since Cholla wasn't compiled with Dual Energy Formalism!\n");
 
#else
  int id, n_cells;
  int neighbor_ids[NDim * 2];
  bool is_real_cell;
 
  // Here we branch and determine some basic details based on the number of dimensions
  // -> there are much more clever ways to do all of this more concisely
  // -> for the uninitiated, the conditions of `constexpr if` statements (within templates) are evaluated
  //    at compile time. In other words, there is no branching
  if constexpr (NDim == 1) {
    n_cells = grid_shape.x();
    id      = threadIdx.x + blockIdx.x * blockDim.x;
 
    int nx = grid_shape.x();
 
    int xid = id;
 
    int imo         = max(xid - 1, n_ghost);
    neighbor_ids[0] = imo;
    int ipo         = min(xid + 1, nx - n_ghost - 1);
    neighbor_ids[1] = ipo;
 
    is_real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost);
 
  } else if constexpr (NDim == 2) {
    n_cells = grid_shape.x() * grid_shape.y();
 
    // this difference in logic is left over from the distant past. when we
    // we needed to use 2D grids to run really large 2D problems, because of
    // CUDA limits on the 1D grid dimensions. This is reflected the name
    // `dim2dgrid` for the grid dimensions in e.g. the VL_2D_cuda.cu functions,
    // versus 'dim1dgrid' for the 3D version, even though then we set up a a 1D
    // grid. In principle as long as we continue to only use 1D grids,
    // Evan thinks that replacing this with blockIdx.x would be fine, since
    // blockIdx.y will always be 0 for a 1D grid,
    int blockId = blockIdx.x + blockIdx.y * gridDim.x;
    id          = threadIdx.x + blockId * blockDim.x;
 
    int nx = grid_shape.x();
    int ny = grid_shape.y();
 
    int yid = id / nx;
    int xid = id - yid * nx;
 
    int imo         = max(xid - 1, n_ghost);
    neighbor_ids[0] = imo + yid * nx;
    int ipo         = min(xid + 1, nx - n_ghost - 1);
    neighbor_ids[1] = ipo + yid * nx;
    int jmo         = max(yid - 1, n_ghost);
    neighbor_ids[2] = xid + jmo * nx;
    int jpo         = min(yid + 1, ny - n_ghost - 1);
    neighbor_ids[3] = xid + jpo * nx;
 
    is_real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost);
 
  } else {  // NDim == 3
    n_cells = grid_shape.x() * grid_shape.y() * grid_shape.z();
    id      = threadIdx.x + blockIdx.x * blockDim.x;
 
    int nx = grid_shape.x();
    int ny = grid_shape.y();
    int nz = grid_shape.z();
 
    int zid = id / (nx * ny);
    int yid = (id - zid * nx * ny) / nx;
    int xid = id - zid * nx * ny - yid * nx;
 
    int imo         = max(xid - 1, n_ghost);
    neighbor_ids[0] = imo + yid * nx + zid * nx * ny;
    int ipo         = min(xid + 1, nx - n_ghost - 1);
    neighbor_ids[1] = ipo + yid * nx + zid * nx * ny;
    int jmo         = max(yid - 1, n_ghost);
    neighbor_ids[2] = xid + jmo * nx + zid * nx * ny;
    int jpo         = min(yid + 1, ny - n_ghost - 1);
    neighbor_ids[3] = xid + jpo * nx + zid * nx * ny;
    int kmo         = max(zid - 1, n_ghost);
    neighbor_ids[4] = xid + yid * nx + kmo * nx * ny;
    int kpo         = min(zid + 1, nz - n_ghost - 1);
    neighbor_ids[5] = xid + yid * nx + kpo * nx * ny;
 
    is_real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost &&
                    zid > n_ghost - 1 && zid < nz - n_ghost);
  }
 
  const Real eta_1 = DE_ETA_1;
  const Real eta_2 = DE_ETA_2;
 
  // threads corresponding to real cells do the calculation
  if (is_real_cell) {
    // every thread collects the conserved variables it needs from global memory
    Real d          = dev_conserved[id];
    Real d_inv      = 1.0 / d;
    Real vx         = dev_conserved[1 * n_cells + id] * d_inv;
    Real vy         = dev_conserved[2 * n_cells + id] * d_inv;
    Real vz         = dev_conserved[3 * n_cells + id] * d_inv;
    Real E          = dev_conserved[4 * n_cells + id];
    Real U_advected = dev_conserved[(n_fields - 1) * n_cells + id];
    Real U_total    = E - 0.5 * d * (vx * vx + vy * vy + vz * vz);
 
    // We will deal with this crashed cell later in a different kernel... (at the time of writing, a 1D
    // simulation always uses floors for this purpose)
    if (Cell_Is_Crashed(d, E)) return;
 
    // find the max nearby total energy (from the local cell and any uncrashed neighbors)
    // -> we take the stance that "crashed" neighbors are unreliable, even if total energy looks ok
    // -> we're making use of a "range-based for loop"
    Real Emax = E;
    for (int neighbor_id : neighbor_ids) {
      Real neighbor_d = dev_conserved[grid_enum::density * n_cells + neighbor_id];
      Real neighbor_E = dev_conserved[grid_enum::Energy * n_cells + neighbor_id];
      Emax            = fmax(Emax, Cell_Is_Crashed(neighbor_d, neighbor_E) ? Emax : neighbor_E);
    }
 
    // Ordinarily, 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)
    // We also explicitly use the "advected" internal energy if the total energy is positive but doesn't
    // exceed kinetic energy (i.e. U_total <= 0). This scenario comes up in simulations with particle-based
    // feedback.
    bool prefer_U_total = (U_total > E * eta_1) or (U_total > Emax * eta_2);
    Real U              = (prefer_U_total and (U_total > 0)) ? U_total : U_advected;
 
    // Optional: Avoid Negative Internal  Energies
    U = fmax(U, (Real)TINY_NUMBER);
 
    // Write Selected internal energy to the GasEnergy array ONLY
    // to avoid mixing updated and non-updated values of E
    // since the Dual Energy condition depends on the neighbour cells
    dev_conserved[(n_fields - 1) * n_cells + id] = U;
  }
#endif /* DE */
}
 
template <int NDim>
__global__ void Sync_Energies(Real *dev_conserved, hydro_utilities::VectorXYZ<int> grid_shape, int n_ghost,
                              int n_fields)
{
  static_assert((NDim == 1) || (NDim == 2) || (NDim == 3), "NDim must be 1, 2, or 3");
 
#ifndef DE
  printf("WARNING: this function isn't usable since Cholla wasn't compiled with Dual Energy Formalism!\n");
 
#else
  int id, n_cells;
  bool is_real_cell;
 
  // Here we branch and determine some basic details based on the number of dimensions
  // -> for the uninitiated, the conditions of `constexpr if` statements (within templates) are evaluated
  //    at compile time. In other words, there is no branching
  if constexpr (NDim == 1) {
    n_cells = grid_shape.x();
 
    // get a global thread ID
    id           = threadIdx.x + blockIdx.x * blockDim.x;
    int xid      = id;
    is_real_cell = (xid > n_ghost - 1 && xid < grid_shape.x() - n_ghost);
 
  } else if constexpr (NDim == 2) {
    n_cells = grid_shape.x() * grid_shape.y();
 
    int nx = grid_shape.x();
    int ny = grid_shape.y();
 
    // get a global thread ID
    int blockId  = blockIdx.x + blockIdx.y * gridDim.x;
    id           = threadIdx.x + blockId * blockDim.x;
    int yid      = id / nx;
    int xid      = id - yid * nx;
    is_real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost);
 
  } else {  // NDim == 3
    n_cells = grid_shape.x() * grid_shape.y() * grid_shape.z();
 
    int nx = grid_shape.x();
    int ny = grid_shape.y();
    int nz = grid_shape.z();
 
    // get a global thread ID
    id           = threadIdx.x + blockIdx.x * blockDim.x;
    int zid      = id / (nx * ny);
    int yid      = (id - zid * nx * ny) / nx;
    int xid      = id - zid * nx * ny - yid * nx;
    is_real_cell = (xid > n_ghost - 1 && xid < nx - n_ghost && yid > n_ghost - 1 && yid < ny - n_ghost &&
                    zid > n_ghost - 1 && zid < nz - n_ghost);
  }
 
  // threads corresponding to real cells do the calculation
  if (is_real_cell) {
    // every thread collects the conserved variables it needs from global memory
    Real d     = dev_conserved[grid_enum::density * n_cells + id];
    Real d_inv = 1.0 / d;
    Real vx    = dev_conserved[grid_enum::momentum_x * n_cells + id] * d_inv;
    Real vy    = dev_conserved[grid_enum::momentum_y * n_cells + id] * d_inv;
    Real vz    = dev_conserved[grid_enum::momentum_z * n_cells + id] * d_inv;
    Real U     = dev_conserved[(n_fields - 1) * n_cells + id];
 
    // Use the previously selected Internal Energy to update the total energy
    dev_conserved[grid_enum::Energy * n_cells + id] = 0.5 * d * (vx * vx + vy * vy + vz * vz) + U;
  }
 
#endif /* DE */
}
 
}  // namespace dual_energy
 
#endif /* DUAL_ENERGY_H */