HydroForce.cpph

Go to the documentation of this file.

#include "equation_of_state.h"
#include "hydro_dimensions.h"
#include "kernel_hydro.h"
 
 
template <typename Particle>
bool swift2::kernels::legacy::forceKernelUpdateParticlePredicate(
  const peano4::datamanagement::CellMarker& marker, const bool& localParticleIsContainedInCell, const Particle& localParticle
) {
  return ::swift2::kernels::localParticleCanBeUpdatedInCellKernel(marker, localParticleIsContainedInCell, localParticle);
}
 
 
template <typename Particle>
bool swift2::kernels::legacy::forceKernelEvaluatePairInteractionPredicate(
  const peano4::datamanagement::CellMarker& marker, const Particle& localParticle, const Particle& activeParticle
) {
  assertion3(marker.isContained(
      localParticle.getX(),
      toolbox::particles::internal::relativeGrabOwnershipSpatialSortingTolerance(marker) * tarch::la::NUMERICAL_ZERO_DIFFERENCE
    ),
    marker.toString(), localParticle.toString(), activeParticle.toString()
  );
  return &localParticle != &activeParticle;
}
 
 
template <typename Particle>
void swift2::kernels::legacy::forceKernel(
  const peano4::datamanagement::CellMarker& marker,
  Particle&                                 localParticle,
  const Particle&                           activeParticle
) {
  assertion2(
    ::swift2::kernels::localParticleCanBeUpdatedInCellKernelFromAnyOtherParticle(
      marker,
      localParticle,
      activeParticle
    ), localParticle.toString(), activeParticle.toString()
  );
  force_kernel(&localParticle, &activeParticle);
}
 
 
template <typename Particle>
void swift2::kernels::legacy::hydro_prepare_force(Particle* localParticle) {
 
  static tarch::logging::Log _log("swift2::kernels::legacy");
 
  // Grab static parameters
  const int    hydroDimensions = localParticle->getHydroDimensions();
  const double alphaAV         = localParticle->getAlphaAV();
  double       hydro_h_max     = std::min(
    localParticle->getSmlMax(),
    localParticle->getSearchRadius() / kernel_gamma
  );
 
  // Retrieve basic properties of local and active particles
  const double rho_i = localParticle->getDensity();
  const double m_i   = localParticle->getMass();
  const double h_i   = localParticle->getSmoothingLength();
  const double u_i   = localParticle->getU();
 
  const double h_inv = 1.0 / h_i;
 
  // This is needed for cosmological integration.
  const double fac_balsara_eps = 1.;
 
  // Compute the norm of the curl
#if Dimensions < 3
  const double curl_v = std::sqrt(
    localParticle->getRot_v() * localParticle->getRot_v()
  );
#else
  const double curl_v = tarch::la::norm2(localParticle->getRot_v());
#endif
 
  /* Compute the norm of div v */
  const double abs_div_v = tarch::la::abs(localParticle->getDiv_v());
 
  /* Compute the pressure */
  const double P_i = eos::gas_pressure_from_internal_energy(rho_i, u_i);
  assertion1(std::isfinite(P_i), localParticle->toString());
  assertion1(P_i > 0., localParticle->toString());
 
  /* Compute the sound speed */
  const double soundspeed = legacy::eos::gas_soundspeed_from_pressure(
    rho_i,
    P_i
  );
  assertion1(std::isfinite(soundspeed), localParticle->toString());
 
  /* Compute the "grad h" term  - Note here that we have \tilde{x}
   * as 1 as we use the local number density to find neighbours. This
   * introduces a j-component that is considered in the force loop,
   * meaning that this cached grad_h_term gives:
   *
   * f_ij = 1.f - grad_h_term_i / m_j */
  const double common_factor = h_i / hydroDimensions
                               / localParticle->getWcount();
 
  float grad_h_term;
  // Ignore changing-kernel effects when h ~= h_max
  if (h_i > 0.9999 * hydro_h_max) {
    grad_h_term = 0.;
    logWarning(
      "hydro_prepare_force",
      "h ~ h_max for particle x=" << localParticle->getX() <<
#if PeanoDebug > 0
        ", partID=" << localParticle->getPartid() <<
#endif
        ", h = " << h_i << ", h_max = " << hydro_h_max
    );
  } else {
    const double grad_W_term = common_factor * localParticle->getWcount_dh();
    if (grad_W_term < -0.9999) {
      /* if we get here, we either had very small neighbour contributions
         (which should be treated as a no neighbour case in the ghost) or
         a very weird particle distribution (e.g. particles sitting on
         top of each other). Either way, we cannot use the normal
         expression, since that would lead to overflow or excessive round
         off and cause excessively high accelerations in the force loop */
      grad_h_term = 0.;
      logWarning(
        "hydro_prepare_force",
        "grad_W_term very small for particle x="
          << localParticle->getX() <<
#if PeanoDebug > 0
          ", partID=" << localParticle->getPartid() <<
#endif
          ", h:" << h_i << ", wcount:" << localParticle->getWcount()
          << ", wcount_dh:" << localParticle->getWcount_dh() << ", grad_W_term:"
          << grad_W_term << ", common_factor:" << common_factor
      );
    } else {
      grad_h_term = common_factor * localParticle->getRho_dh()
                    / (1. + grad_W_term);
    }
  }
 
 
  /* Compute the Balsara Switch */
  const double balsara
    = alphaAV * abs_div_v
      / (abs_div_v + curl_v + 0.0001 * fac_balsara_eps * soundspeed * h_inv);
 
  assertion6(
    std::isfinite(balsara),
    localParticle->toString(),
    balsara,
    abs_div_v,
    curl_v,
    soundspeed,
    alphaAV
  );
 
  /* Update variables */
  localParticle->setF(grad_h_term);
  localParticle->setPressure(P_i);
  localParticle->setSoundSpeed(soundspeed);
  localParticle->setBalsara(balsara);
 
#if PeanoDebug > 0
  localParticle->setForceNeighbourCount(0);
#endif
}
 
 
template <typename Particle>
void swift2::kernels::legacy::hydro_reset_acceleration(Particle* localParticle
) {
 
  /* Reset the acceleration*/
  localParticle->setA(0.0);
 
  /* Reset the time derivatives */
  localParticle->setUDot(0.0);
  localParticle->setHDot(0.0);
  localParticle->setV_sig_AV(2.0 * localParticle->getSoundSpeed());
}
 
 
template <typename Particle>
void swift2::kernels::legacy::hydro_end_force(Particle* localParticle) {
 
  static tarch::logging::Log _log("swift2::kernels::legacy");
 
  const int    hydroDimensions    = localParticle->getHydroDimensions();
  const double hydroDimensionsInv = 1. / (double)hydroDimensions;
 
  localParticle->setHDot(
    localParticle->getHDot() * localParticle->getSmoothingLength()
    * hydroDimensionsInv
  );
  assertion2(
    std::isfinite(localParticle->getHDot()),
    localParticle->toString(),
    localParticle->getSmoothingLength()
  );
 
#if PeanoDebug > 0
  // During the smoothing length computation, we complete a neighbour loop
  // while assuming we have the correct sml. In the final step, we compute
  // the new sml and compare it to the one we used. If the difference is
  // small enough, as set by the convergence criterion, we consider the
  // sml converged. However, it is imporant to note that we use "small
  // enough", not "zero". So the sml may be slightly different after the
  // sml computation iteration completes.
  // Down the line, this can lead to scenarios where we interacted with
  // more particles during the density loop than the force loop. Provided
  // the convergence criterion has been set adequately small, the difference
  // in results should be negligible. However, since for a fixed smoothing
  // length one expects equal or more neighbour interactions during a force
  // loop since we take into account either a particle's own radius or the
  // radius of the interacting neighbour, as opposed to only the particle's
  // own radius during the density loop, this expectation is no longer strictly
  // satisfied.
  // Previously, I used to have an assertion here to make sure the interaction
  // counts satisfy that expectation. For the reason explained above, that
  // doesn't work. So instead, we just check for it and raise a warning, since
  // it *could* be an indication of something going wrong. Typically, one can
  // expect a 1-2 particle difference. If it's larger than that, there is
  // most likely an issue present.
  if (localParticle->getDensityNeighbourCount() > localParticle->getForceNeighbourCount()) {
    logWarning(
      "hydro_end_force()",
      "Particle interacted with fewer particles in Force than in Density loop"
        << localParticle->toString()
    )
  }
#endif
}
 
 
template <typename Particle>
void swift2::kernels::legacy::force_kernel(
  Particle*       localParticle,
  const Particle* activeParticle
) {
 
  // macro expansions from swift2/legacy/kernels are cumbersome otherwise.
  using namespace kernelHydro;
 
  // Grab static parameters
  const int    hydroDimensions = localParticle->getHydroDimensions();
  const double alphaAV         = localParticle->getAlphaAV();
  const double betaAV          = localParticle->getBetaAV();
 
  // No cosmology yet: a^{(3 gamma - 5 )/ 2} = 1.
  const double fac_mu    = 1.;
  const double a2_Hubble = 0.; // a**2 * H, where H = da/dt / a. = 0 without
                               // cosmo
 
  // Distance between the two particles
  const tarch::la::Vector<Dimensions, double>
    dx = localParticle->getX() - activeParticle->getX();
 
  const double r     = tarch::la::norm2(dx);
  const double r_inv = r ? 1.0 / r : 0.0;
 
  // Get actual interaction radius
  const double h_i    = localParticle->getSmoothingLength();
  const double iactRi = kernel_gamma * h_i;
  const double h_j    = activeParticle->getSmoothingLength();
  const double iactRj = kernel_gamma * h_j;
 
  if ((tarch::la::smaller(r, iactRi) or (tarch::la::smaller(r, iactRj))) and tarch::la::greater(r, 0.0)) {
 
#if PeanoDebug > 0
    localParticle->setForceNeighbourCount(
      localParticle->getForceNeighbourCount() + 1
    );
#endif
 
    // Retrieve basic properties of local and active particles
    const double m_i = localParticle->getMass();
    const double m_j = activeParticle->getMass();
 
    const double rho_i = localParticle->getDensity();
    const double rho_j = activeParticle->getDensity();
 
    const double P_i = localParticle->getPressure();
    const double P_j = activeParticle->getPressure();
 
    // Get kernel for hi
    const double h_i_inv     = 1. / h_i;
    const double h_i_inv_dp1 = swift2::kernels::legacy::hydroDimensions::
      pow_dimension_plus_one(h_i_inv); // 1/h^{dim + 1}
    const double qi = r * h_i_inv;
    double       wi, dwi_dqi;
    // Evaluate kernel and its derivative for particle a
    kernel_deval(qi, wi, dwi_dqi);
    const double dwi_dr = dwi_dqi * h_i_inv_dp1;
 
    // Get kernel for hj
    const double h_j_inv     = 1. / h_j;
    const double h_j_inv_dp1 = swift2::kernels::legacy::hydroDimensions::
      pow_dimension_plus_one(h_j_inv); // 1/h^{dim + 1}
    const double qj = r * h_j_inv;
    double       wj, dwj_dqj;
    // Evaluate kernel and its derivative for particle b
    kernel_deval(qj, wj, dwj_dqj);
    const double dwj_dr = dwj_dqj * h_j_inv_dp1;
 
    // Variable smoothing length term
 
    // SWIFT "f"
    const double f_ij = 1.0 - localParticle->getF() / m_j;
    const double f_ji = 1.0 - activeParticle->getF() / m_i;
 
    // Compute gradient terms
    const double P_over_rho2_i = P_i / (rho_i * rho_i) * f_ij;
    const double P_over_rho2_j = P_j / (rho_j * rho_j) * f_ji;
 
    // Compute dv dot r.
    double dvdr = tarch::la::dot(
      (localParticle->getV() - activeParticle->getV()),
      dx
    );
 
    // Add Hubble flow: Skipped for now
    const double dvdr_Hubble = dvdr + a2_Hubble * r * r;
 
    /* Are the particles moving towards each others ? */
    const double omega_ij = std::min(dvdr_Hubble, 0.);
    const double mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
    /* Compute sound speeds and signal velocity */
    const double v_sig = localParticle->getSoundSpeed()
                         + activeParticle->getSoundSpeed() - betaAV * mu_ij;
 
    /* Balsara term */
    const double balsara_i = localParticle->getBalsara();
    const double balsara_j = activeParticle->getBalsara();
 
    /* Construct the full viscosity term */
    const double rho_ij = 0.5 * (rho_i + rho_j);
    const double visc = -0.25 * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
    /* Convolve with the kernel */
    const double visc_acc_term = 0.5 * visc * (dwi_dr * f_ij + dwj_dr * f_ji)
                                 * r_inv;
 
    /* SPH acceleration term */
    const double sph_acc_term = (P_over_rho2_i * dwi_dr + P_over_rho2_j * dwj_dr)
                                * r_inv;
 
    // Assemble the acceleration
    const double acc = sph_acc_term + visc_acc_term;
 
    assertion9(
      std::isfinite(visc_acc_term),
      visc_acc_term,
      localParticle->toString(),
      activeParticle->toString(),
      dwi_dr,
      f_ij,
      dwj_dr,
      f_ji,
      r_inv,
      visc
    );
    assertion3(
      std::isfinite(sph_acc_term),
      sph_acc_term,
      localParticle->toString(),
      activeParticle->toString()
    );
 
    // Use the force
    // NOTE: This is a_hydro in SWIFT
    localParticle->setA(localParticle->getA() - m_j * acc * dx);
 
    /* Get the time derivative for u. */
    const double sph_du_term_i = P_over_rho2_i * dvdr * r_inv * dwi_dr;
 
    /* Viscosity term */
    const double visc_du_term = 0.5 * visc_acc_term * dvdr_Hubble;
 
    // Assemble the energy equation term
    const double du_dt_i = sph_du_term_i + visc_du_term;
 
    // Internal energy time derivative
    localParticle->setUDot(localParticle->getUDot() + m_j * du_dt_i);
 
    /* Get the time derivative for h. */
    localParticle->setHDot(
      localParticle->getHDot() - m_j * dvdr * r_inv / rho_j * dwi_dr * f_ij
    );
 
    // Update the signal velocity
    localParticle->setV_sig_AV(std::max(localParticle->getV_sig_AV(), v_sig));
 
    assertion6(
      not std::isnan(localParticle->getHDot()),
      localParticle->toString(),
      dvdr,
      r_inv,
      dwi_dr,
      f_ij,
      f_ji
    );
    assertion6(
      not tarch::la::isEntryNan(localParticle->getA()),
      localParticle->toString(),
      sph_acc_term,
      visc_acc_term,
      f_ij,
      f_ji,
      activeParticle->toString()
    );
  }
}