FluidDynamics.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FluidDynamics.h

/*******************************************************************************
 * Copyright (c) The JETSCAPE Collaboration, 2018
 *
 * Modular, task-based framework for simulating all aspects of heavy-ion collisions
 *
 * For the list of contributors see AUTHORS.
 *
 * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
 *
 * or via email to bugs.jetscape@gmail.com
 *
 * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
 * See COPYING for details.
 ******************************************************************************/
// This is a general basic class for hydrodynamics

#ifndef FLUIDDYNAMICS_H
#define FLUIDDYNAMICS_H

#include <memory>
#include <vector>
#include <array>
#include <cstring>
#include <stdexcept>
#include <cmath>
#include <iostream>
#include <map>

#include "InitialState.h"
#include "JetScapeModuleBase.h"
#include "PreequilibriumDynamics.h"
#include "RealType.h"
#include "FluidCellInfo.h"
#include "FluidEvolutionHistory.h"
#include "LiquefierBase.h"
#include "SurfaceCellInfo.h"

namespace Jetscape {

enum HydroStatus { NOT_START, INITIALIZED, EVOLVING, FINISHED, ERROR };

class Parameter {
public:
  char *hydro_input_filename;
};

class FluidDynamics : public JetScapeModuleBase {
protected:
  // record hydro start and end proper time [fm/c]
  Jetscape::real hydro_tau_0, hydro_tau_max;
  // record hydro freeze out temperature [GeV]
  Jetscape::real hydro_freeze_out_temperature;
  // record hydro running status
  HydroStatus hydro_status;

  // add initial state shared pointer
  std::shared_ptr<InitialState> ini;
  std::shared_ptr<PreequilibriumDynamics> pre_eq_ptr;

  double eta;
  bool boost_invariant_;
  Parameter parameter_list;

  // How to store this data? In memory or hard disk?
  // 3D hydro may eat out the memory,
  // for large dataset, std::deque is better than std::vector.
  EvolutionHistory bulk_info;

  std::weak_ptr<LiquefierBase> liquefier_ptr;

public:
  FluidDynamics();

  virtual ~FluidDynamics();

  virtual void Init();

  virtual void Exec();

  virtual void Clear();

  Parameter &GetParameterList() { return (parameter_list); }

  virtual void CollectHeader(weak_ptr<JetScapeWriter> w);

  virtual double GetEventPlaneAngle() { return (-999); }

  virtual void GetTemperature(double t, double x, double y, double z,
                              double &mT) {
    mT = GetTemperature(t, x, y, z);
  }

  virtual void GetHydroCell(double t, double x, double y, double z,
                            std::unique_ptr<FluidCellInfo> &fCell) {
    GetHydroInfo(t, x, y, z, fCell);
  }

  // currently we have no standard for passing configurations
  // pure virtual function; to be implemented by users
  // should make it easy to save evolution history to bulk_info
  virtual void InitializeHydro(Parameter parameter_list){};

  virtual void EvolveHydro(){};

  int GetHydroStatus() const { return (hydro_status); }

  void StoreHydroEvolutionHistory(
      std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
    bulk_info.data.push_back(*fluid_cell_info_ptr);
  }

  void clear_up_evolution_data() { bulk_info.clear_up_evolution_data(); }

  void GetHydroStartTime(double &tau0) { tau0 = hydro_tau_0; }

  Jetscape::real GetHydroEndTime() const { return (hydro_tau_max); }
  Jetscape::real GetHydroFreezeOutTemperature() const {
    return (hydro_freeze_out_temperature);
  }

  virtual void
  GetHydroInfo(Jetscape::real t, Jetscape::real x, Jetscape::real y,
               Jetscape::real z,
               std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
    if (hydro_status != FINISHED || bulk_info.data.size() == 0) {
      throw std::runtime_error("Hydro evolution is not finished "
                               "or EvolutionHistory is empty");
    }
    // judge whether to use 2D interpolation or 3D interpolation
    if (!bulk_info.tau_eta_is_tz) {
      Jetscape::real tau = std::sqrt(t * t - z * z);
      Jetscape::real eta = 0.5 * (std::log(t + z) - std::log(t - z));
      bulk_info.CheckInRange(tau, x, y, eta);
      //return bulk_info.get(tau, x, y, eta);
    } else {
      bulk_info.CheckInRange(t, x, y, z);
      //return bulk_info.get(t, x, y, z);
    }
  }

  // this function print out the information of the fluid cell to the screen
  void PrintFluidCellInformation(FluidCellInfo *fluid_cell_info_ptr);

  // this function returns hypersurface for Cooper-Frye or recombination
  // the detailed implementation is left to the hydro developper
  void FindAConstantTemperatureSurface(
          Jetscape::real T_sw, std::vector<SurfaceCellInfo> &surface_cells);

  // all the following functions will call function GetHydroInfo()
  // to get thermaldynamic and dynamical information at a space-time point
  // (time, x, y, z)

  virtual Jetscape::real GetEnergyDensity(Jetscape::real time, Jetscape::real x,
                                          Jetscape::real y, Jetscape::real z);

  virtual Jetscape::real GetEntropyDensity(Jetscape::real time,
                                           Jetscape::real x, Jetscape::real y,
                                           Jetscape::real z);

  virtual Jetscape::real GetTemperature(Jetscape::real time, Jetscape::real x,
                                        Jetscape::real y, Jetscape::real z);

  virtual Jetscape::real GetQgpFraction(Jetscape::real time, Jetscape::real x,
                                        Jetscape::real y, Jetscape::real z);

  // These have no default implementation
  // /** @return 3-component (vx,vy,vz) fluid velocity at point (t or tau, x, y, z or eta).
  //     @param time Time or tau coordinate.
  //     @param x Space coordinate.
  //     @param y Space coordinate.
  //     @param z Space or eta coordinate.
  // */
  // virtual real3 Get3FluidVelocity(Jetscape::real time, Jetscape::real x, Jetscape::real y, Jetscape::real z)=0;

  // /** @return 4-component fluid velocity at point (t or tau, x, y, zor eta).
  //     @param time Time or tau coordinate.
  //     @param x Space coordinate.
  //     @param y Space coordinate.
  //     @param z Space or eta coordinate.
  // */
  // virtual real4 Get4FluidVelocity(Jetscape::real time, Jetscape::real x, Jetscape::real y, Jetscape::real z);

  // /** @return Net baryon density at point (t or tau, x, y, z or eta).
  //     @param time Time or tau coordinate.
  //     @param x Space coordinate.
  //     @param y Space coordinate.
  //     @param z Space or eta coordinate.
  // */
  // virtual Jetscape::real GetNetBaryonDensity(Jetscape::real time, Jetscape::real x, Jetscape::real y, Jetscape::real z);

  // /** @return Net charge density at point (t or tau, x, y, z or eta).
  //    @param time Time or tau coordinate.
  //    @param x Space coordinate.
  //    @param y Space coordinate.
  //    @param z Space or eta coordinate.
  // */
  // virtual Jetscape::real GetNetChargeDensity(Jetscape::real time, Jetscape::real x, Jetscape::real y, Jetscape::real z);

  virtual void add_a_liquefier(std::shared_ptr<LiquefierBase> new_liquefier) {
    liquefier_ptr = new_liquefier;
  }

  void get_source_term(Jetscape::real tau, Jetscape::real x, Jetscape::real y,
                       Jetscape::real eta,
                       std::array<Jetscape::real, 4> jmu) const;

  virtual void UpdateEnergyDeposit(int t, double edop);
  virtual void GetEnergyDensity(int t, double &edensity) { edensity = 0.0; }

  // get a reference to the bulk_info object
  const EvolutionHistory& get_bulk_info() const { return bulk_info; }

}; // end class FluidDynamics

} // end namespace Jetscape

#endif // FLUIDDYNAMICS_H