density.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2022,2024
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
#ifndef SRC_INCLUDE_SMASH_DENSITY_H_
#define SRC_INCLUDE_SMASH_DENSITY_H_
 
#include <iostream>
#include <tuple>
#include <typeinfo>
#include <utility>
#include <vector>
 
#include "energymomentumtensor.h"
#include "experimentparameters.h"
#include "forwarddeclarations.h"
#include "fourvector.h"
#include "lattice.h"
#include "particledata.h"
#include "particles.h"
#include "pdgcode.h"
#include "threevector.h"
 
namespace smash {
static constexpr int LDensity = LogArea::Density::id;
 
std::ostream &operator<<(std::ostream &os, DensityType dt);
 
double density_factor(const ParticleType &type, DensityType dens_type);
 
inline double smearing_factor_norm(const double two_sigma_sqr) {
  const double tmp = two_sigma_sqr * M_PI;
  return tmp * std::sqrt(tmp);
}
 
inline double smearing_factor_rcut_correction(const double rcut_in_sigma) {
  const double x = rcut_in_sigma / std::sqrt(2.0);
  return -2.0 / std::sqrt(M_PI) * x * std::exp(-x * x) + std::erf(x);
}
 
class DensityParameters {
 public:
  DensityParameters(const ExperimentParameters &par)  // NOLINT
      : sig_(par.gaussian_sigma),
        r_cut_(par.gauss_cutoff_in_sigma * par.gaussian_sigma),
        ntest_(par.testparticles),
        nensembles_(par.n_ensembles),
        derivatives_(par.derivatives_mode),
        rho_derivatives_(par.rho_derivatives_mode),
        smearing_(par.smearing_mode),
        central_weight_(par.discrete_weight),
        triangular_range_(par.triangular_range),
        only_participants_(par.only_participants) {
    r_cut_sqr_ = r_cut_ * r_cut_;
    const double two_sig_sqr = 2 * sig_ * sig_;
    two_sig_sqr_inv_ = 1. / two_sig_sqr;
    const double norm = smearing_factor_norm(two_sig_sqr);
    const double corr_factor =
        smearing_factor_rcut_correction(par.gauss_cutoff_in_sigma);
    norm_factor_sf_ = 1. / (norm * ntest_ * nensembles_ * corr_factor);
  }
  int ntest() const { return ntest_; }
  int nensembles() const { return nensembles_; }
  DerivativesMode derivatives() const { return derivatives_; }
  RestFrameDensityDerivativesMode rho_derivatives() const {
    return rho_derivatives_;
  }
  SmearingMode smearing() const { return smearing_; }
  double central_weight() const { return central_weight_; }
  double triangular_range() const { return triangular_range_; }
  double r_cut() const { return r_cut_; }
  double r_cut_sqr() const { return r_cut_sqr_; }
  double two_sig_sqr_inv() const { return two_sig_sqr_inv_; }
  double norm_factor_sf() const { return norm_factor_sf_; }
  bool only_participants() const { return only_participants_; }
 
 private:
  const double sig_;
  const double r_cut_;
  double r_cut_sqr_;
  double two_sig_sqr_inv_;
  double norm_factor_sf_;
  const int ntest_;
  const int nensembles_;
  const DerivativesMode derivatives_;
  const RestFrameDensityDerivativesMode rho_derivatives_;
  const SmearingMode smearing_;
  const double central_weight_;
  const double triangular_range_;
  bool only_participants_;
};
 
std::pair<double, ThreeVector> unnormalized_smearing_factor(
    const ThreeVector &r, const FourVector &p, const double m_inv,
    const DensityParameters &dens_par, const bool compute_gradient = false);
 
std::tuple<double, FourVector, ThreeVector, ThreeVector, FourVector, FourVector,
           FourVector, FourVector>
current_eckart(const ThreeVector &r, const ParticleList &plist,
               const DensityParameters &par, DensityType dens_type,
               bool compute_gradient, bool smearing);
std::tuple<double, FourVector, ThreeVector, ThreeVector, FourVector, FourVector,
           FourVector, FourVector>
current_eckart(const ThreeVector &r, const Particles &plist,
               const DensityParameters &par, DensityType dens_type,
               bool compute_gradient, bool smearing);
 
class DensityOnLattice {
 public:
  DensityOnLattice()
      : jmu_pos_(FourVector()),
        jmu_neg_(FourVector()),
        djmu_dxnu_({FourVector(), FourVector(), FourVector(), FourVector()}),
        drho_dxnu_(FourVector()) {}
 
  void add_particle(const ParticleData &part, double FactorTimesSf) {
    const FourVector part_four_velocity = FourVector(1.0, part.velocity());
    if (FactorTimesSf > 0.0) {
      jmu_pos_ += part_four_velocity * FactorTimesSf;
    } else {
      jmu_neg_ += part_four_velocity * FactorTimesSf;
    }
  }
 
  void add_particle_for_derivatives(const ParticleData &part, double factor,
                                    ThreeVector sf_grad) {
    const FourVector PartFourVelocity = FourVector(1.0, part.velocity());
    for (int k = 1; k <= 3; k++) {
      djmu_dxnu_[k] += factor * PartFourVelocity * sf_grad[k - 1];
      djmu_dxnu_[0] -=
          factor * PartFourVelocity * sf_grad[k - 1] * part.velocity()[k - 1];
    }
  }
 
  double rho(const double norm_factor = 1.0) {
    return (jmu_pos_.abs() - jmu_neg_.abs()) * norm_factor;
  }
 
  ThreeVector curl_vecj(const double norm_factor = 1.0) {
    ThreeVector curl_vec_j = ThreeVector();
    curl_vec_j.set_x1(djmu_dxnu_[2].x3() - djmu_dxnu_[3].x2());
    curl_vec_j.set_x2(djmu_dxnu_[3].x1() - djmu_dxnu_[1].x3());
    curl_vec_j.set_x3(djmu_dxnu_[1].x2() - djmu_dxnu_[2].x1());
    curl_vec_j *= norm_factor;
    return curl_vec_j;
  }
 
  ThreeVector grad_j0(const double norm_factor = 1.0) {
    ThreeVector j0_grad = ThreeVector();
    for (int i = 1; i < 4; i++) {
      j0_grad[i - 1] = djmu_dxnu_[i].x0() * norm_factor;
    }
    return j0_grad;
  }
 
  ThreeVector dvecj_dt(const double norm_factor = 1.0) {
    return djmu_dxnu_[0].threevec() * norm_factor;
  }
 
  FourVector jmu_net() const { return jmu_pos_ + jmu_neg_; }
 
  void add_to_jmu_pos(FourVector additional_jmu_B) {
    jmu_pos_ += additional_jmu_B;
  }
 
  void add_to_jmu_neg(FourVector additional_jmu_B) {
    jmu_neg_ += additional_jmu_B;
  }
 
  FourVector drho_dxnu() const { return drho_dxnu_; }
 
  std::array<FourVector, 4> djmu_dxnu() const { return djmu_dxnu_; }
 
  ThreeVector grad_rho_cross_vecj() const {
    const ThreeVector grad_rho = drho_dxnu_.threevec();
    const ThreeVector vecj = jmu_net().threevec();
    const ThreeVector Drho_cross_vecj = grad_rho.cross_product(vecj);
 
    return Drho_cross_vecj;
  }
 
  void overwrite_djmu_dt_to_zero() {
    djmu_dxnu_[0] = FourVector(0.0, 0.0, 0.0, 0.0);
  }
 
  void overwrite_drho_dt_to_zero() { drho_dxnu_[0] = 0.0; }
 
  void overwrite_drho_dxnu(FourVector computed_drho_dxnu) {
    drho_dxnu_ = computed_drho_dxnu;
  }
 
  void overwrite_djmu_dxnu(FourVector djmu_dt, FourVector djmu_dx,
                           FourVector djmu_dy, FourVector djmu_dz) {
    djmu_dxnu_[0] = djmu_dt;
    djmu_dxnu_[1] = djmu_dx;
    djmu_dxnu_[2] = djmu_dy;
    djmu_dxnu_[3] = djmu_dz;
  }
 
 private:
  FourVector jmu_pos_;
  FourVector jmu_neg_;
  std::array<FourVector, 4> djmu_dxnu_;
  FourVector drho_dxnu_;
};
 
typedef RectangularLattice<DensityOnLattice> DensityLattice;
 
template <typename T>
void update_lattice_with_list_of_particles(RectangularLattice<T> *lat,
                                           const LatticeUpdate update,
                                           const DensityType dens_type,
                                           const DensityParameters &par,
                                           const ParticleList &plist,
                                           const bool compute_gradient,
                                           const bool lattice_reset = true) {
  // Do not proceed if lattice does not exists/update not required
  if (lat == nullptr || lat->when_update() != update) {
    return;
  }
  if (lattice_reset) {
    lat->reset();
  }
  for (const ParticleData &part : plist) {
    if (par.only_participants()) {
      // if this conditions holds, the hadron is a spectator
      if (part.get_history().collisions_per_particle == 0) {
        continue;
      }
    }
    const double dens_factor = density_factor(part.type(), dens_type);
    if (std::abs(dens_factor) < really_small) {
      continue;
    }
    const FourVector p_mu = part.momentum();
    const ThreeVector pos = part.position().threevec();
 
    // act accordingly to which smearing is used
    if (par.smearing() == SmearingMode::CovariantGaussian) {
      // get the normalization factor for the covariant Gaussian smearing
      const double norm_factor_gaus = par.norm_factor_sf();
      const double m = p_mu.abs();
      if (unlikely(m < really_small)) {
        logg[LDensity].warn("Gaussian smearing is undefined for momentum ",
                            p_mu);
        continue;
      }
      const double m_inv = 1.0 / m;
 
      // unweighted contribution to density
      const double common_weight = dens_factor * norm_factor_gaus;
      lat->iterate_in_cube(
          pos, par.r_cut(), [&](T &node, int ix, int iy, int iz) {
            // find the weight for smearing
            const ThreeVector r = lat->cell_center(ix, iy, iz);
            const auto sf = unnormalized_smearing_factor(pos - r, p_mu, m_inv,
                                                         par, compute_gradient);
            node.add_particle(part, sf.first * common_weight);
            if (par.derivatives() == DerivativesMode::CovariantGaussian) {
              node.add_particle_for_derivatives(part, dens_factor,
                                                sf.second * norm_factor_gaus);
            }
          });
    } else if (par.smearing() == SmearingMode::Discrete) {
      // get the volume of the cell and weights for discrete smearing
      const double V_cell = (lat->cell_sizes())[0] * (lat->cell_sizes())[1] *
                            (lat->cell_sizes())[2];
      // weights for coarse smearing
      const double big = par.central_weight();
      const double small = (1.0 - big) / 6.0;
      // unweighted contribution to density
      const double common_weight =
          dens_factor / (par.ntest() * par.nensembles() * V_cell);
      lat->iterate_nearest_neighbors(
          pos, [&](T &node, int iterated_index, int center_index) {
            node.add_particle(
                part, common_weight *
                          // the contribution to density is weighted depending
                          // on what node it is added to
                          (iterated_index == center_index ? big : small));
          });
    } else if (par.smearing() == SmearingMode::Triangular) {
      // get the radii for triangular smearing
      const std::array<double, 3> triangular_radius = {
          par.triangular_range() * (lat->cell_sizes())[0],
          par.triangular_range() * (lat->cell_sizes())[1],
          par.triangular_range() * (lat->cell_sizes())[2]};
      const double prefactor_triangular =
          1.0 /
          (par.ntest() * par.nensembles() * triangular_radius[0] *
           triangular_radius[0] * triangular_radius[1] * triangular_radius[1] *
           triangular_radius[2] * triangular_radius[2]);
      // unweighted contribution to density
      const double common_weight = dens_factor * prefactor_triangular;
      lat->iterate_in_rectangle(
          pos, triangular_radius, [&](T &node, int ix, int iy, int iz) {
            // compute the position of the node
            const ThreeVector cell_center = lat->cell_center(ix, iy, iz);
            // compute smearing weight
            const double weight_x =
                triangular_radius[0] - std::abs(cell_center[0] - pos[0]);
            const double weight_y =
                triangular_radius[1] - std::abs(cell_center[1] - pos[1]);
            const double weight_z =
                triangular_radius[2] - std::abs(cell_center[2] - pos[2]);
            // add the contribution to the node
            node.add_particle(part,
                              common_weight * weight_x * weight_y * weight_z);
          });
    }
  }
}
 
template <typename T>
void update_lattice_accumulating_ensembles(
    RectangularLattice<T> *lat, const LatticeUpdate update,
    const DensityType dens_type, const DensityParameters &par,
    const std::vector<Particles> &ensembles, const bool compute_gradient) {
  // Do not proceed if lattice does not exists/update not required
  if (lat == nullptr || lat->when_update() != update) {
    return;
  }
  lat->reset();
  for (const Particles &particles : ensembles) {
    update_lattice_with_list_of_particles(lat, update, dens_type, par,
                                          particles.copy_to_vector(),
                                          compute_gradient, false);
  }
}
 
void update_lattice(
    RectangularLattice<DensityOnLattice> *lat,
    RectangularLattice<FourVector> *old_jmu,
    RectangularLattice<FourVector> *new_jmu,
    RectangularLattice<std::array<FourVector, 4>> *four_grad_lattice,
    const LatticeUpdate update, const DensityType dens_type,
    const DensityParameters &par, const std::vector<Particles> &ensembles,
    const double time_step, const bool compute_gradient);
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_DENSITY_H_