src/hydro/CLViscWrapper.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #include <stdio.h>
0017 #include <sys/stat.h>
0018 #include <MakeUniqueHelper.h>
0019
0020 #include <string>
0021 #include <vector>
0022
0023 #include "JetScapeLogger.h"
0024 #include "CLViscWrapper.h"
0025
0026 using namespace Jetscape;
0027
0028 // Register the module with the base class
0029 RegisterJetScapeModule<CLVisc> CLVisc::reg("CLVisc");
0030
0031 CLVisc::CLVisc() {
0032   hydro_status = NOT_START;
0033   SetId("CLVisc");
0034 }
0035
0036 CLVisc::~CLVisc() {}
0037
0038 void CLVisc::InitializeHydro(Parameter parameter_list) {
0039   JSINFO << "Initialize CLVisc ...";
0040   VERBOSE(8);
0041
0042   std::string s = GetXMLElementText({"Hydro", "CLVisc", "name"});
0043   JSDEBUG << s << " to be initialized ...";
0044
0045   clvisc::Config cfg;
0046
0047   std::string device_type =
0048       GetXMLElementText({"Hydro", "CLVisc", "device_type"});
0049
0050   if (device_type == "cpu" || device_type == "CPU") {
0051     cfg.block_size = GetXMLElementInt({"Hydro", "CLVisc", "cpu_block_size"});
0052   } else {
0053     cfg.block_size = GetXMLElementInt({"Hydro", "CLVisc", "gpu_block_size"});
0054   }
0055   int device_id = GetXMLElementInt({"Hydro", "CLVisc", "device_id"});
0056   cfg.etaos_xmin = GetXMLElementInt({"Hydro", "CLVisc", "etaos_xmin"});
0057   cfg.etaos_ymin = GetXMLElementInt({"Hydro", "CLVisc", "etaos_ymin"});
0058   cfg.etaos_left_slop =
0059       GetXMLElementInt({"Hydro", "CLVisc", "etaos_left_slop"});
0060   cfg.etaos_right_slop =
0061       GetXMLElementInt({"Hydro", "CLVisc", "etaos_right_slop"});
0062   cfg.result_directory =
0063       GetXMLElementInt({"Hydro", "CLVisc", "result_directory"});
0064   doCooperFrye =
0065       GetXMLElementInt({"Hydro", "CLVisc", "Perform_CooperFrye_Feezeout"});
0066   cfg.tau0 = GetXMLElementDouble({"Hydro", "CLVisc", "tau0"});
0067   cfg.dt = GetXMLElementDouble({"Hydro", "CLVisc", "dtau"});
0068   cfg.ntskip = GetXMLElementDouble({"Hydro", "CLVisc", "ntau_skip"});
0069   cfg.nxskip = GetXMLElementDouble({"Hydro", "CLVisc", "nx_skip"});
0070   cfg.nyskip = GetXMLElementDouble({"Hydro", "CLVisc", "ny_skip"});
0071   cfg.nzskip = GetXMLElementDouble({"Hydro", "CLVisc", "netas_skip"});
0072
0073   cfg.dx = ini->GetXStep();
0074   cfg.dy = ini->GetYStep();
0075   cfg.dz = ini->GetZStep();
0076   cfg.nx = ini->GetXSize();
0077   cfg.ny = ini->GetYSize();
0078   cfg.nz = ini->GetZSize();
0079
0080   hydro_ = std::unique_ptr<clvisc::CLVisc>(
0081       new clvisc::CLVisc(cfg, device_type, device_id));
0082   initial_condition_scale_factor =
0083       GetXMLElementDouble({"Hydro", "CLVisc", "scale_factor"});
0084 }
0085
0086 void CLVisc::EvolveHydro() {
0087   VERBOSE(8);
0088   JSINFO << "Initialize density profiles in CLVisc ...";
0089   std::vector<double> entropy_density = ini->GetEntropyDensityDistribution();
0090   double dx = ini->GetXStep();
0091   if (pre_eq_ptr == nullptr) {
0092     if (initial_condition_scale_factor == 1.0) {
0093       hydro_->read_ini(entropy_density);
0094     } else {
0095       std::for_each(
0096           entropy_density.begin(), entropy_density.end(),
0097           [&](double &sd) { sd = initial_condition_scale_factor * sd; });
0098       hydro_->read_ini(entropy_density);
0099     }
0100   } else {
0101     std::vector<double> vx_, vy_, vz_;
0102     for (size_t idx = 0; idx < pre_eq_ptr->ux_.size(); idx++) {
0103       vx_.push_back(pre_eq_ptr->ux_.at(idx) / pre_eq_ptr->utau_.at(idx));
0104       vy_.push_back(pre_eq_ptr->uy_.at(idx) / pre_eq_ptr->utau_.at(idx));
0105       vz_.push_back(pre_eq_ptr->ueta_.at(idx) / pre_eq_ptr->utau_.at(idx));
0106     }
0107
0108     hydro_->read_ini(pre_eq_ptr->e_, vx_, vy_, vz_, pre_eq_ptr->pi00_,
0109                      pre_eq_ptr->pi01_, pre_eq_ptr->pi02_, pre_eq_ptr->pi03_,
0110                      pre_eq_ptr->pi11_, pre_eq_ptr->pi12_, pre_eq_ptr->pi13_,
0111                      pre_eq_ptr->pi22_, pre_eq_ptr->pi23_, pre_eq_ptr->pi33_);
0112   }
0113
0114   hydro_status = INITIALIZED;
0115   if (hydro_status == INITIALIZED) {
0116     JSINFO << "running CLVisc ...";
0117     hydro_->evolve();
0118     hydro_status = FINISHED;
0119
0120     auto cfg = hydro_->get_config();
0121     float tau_min = cfg.tau0;
0122     float dtau = cfg.dt * cfg.ntskip;
0123     float dx = cfg.dx * cfg.nxskip;
0124     float dy = cfg.dy * cfg.nyskip;
0125     float detas = cfg.dz * cfg.nzskip;
0126     float x_min = -0.5f * cfg.nx * cfg.dx;
0127     float y_min = -0.5f * cfg.ny * cfg.dy;
0128     float etas_min = -0.5f * cfg.nz * cfg.dz;
0129     int nx = int(floor((cfg.nx - 1) / cfg.nxskip)) + 1;
0130     int ny = int(floor((cfg.ny - 1) / cfg.nyskip)) + 1;
0131     int netas = int(floor((cfg.nz - 1) / cfg.nzskip)) + 1;
0132
0133     bulk_info.FromVector(hydro_->bulkinfo_.get_data(),
0134                          hydro_->bulkinfo_.get_data_info(), tau_min, dtau,
0135                          x_min, dx, nx, y_min, dy, ny, etas_min, detas, netas,
0136                          false);
0137     //hydro_->bulkinfo_.save("bulk_data.csv");
0138   }
0139   if (hydro_status == FINISHED && doCooperFrye == 1) {
0140     JSINFO << "Cooper Frye not implemented yet";
0141   }
0142 }
0143
0144 void CLVisc::GetHydroInfo(Jetscape::real t, Jetscape::real x, Jetscape::real y,
0145                           Jetscape::real z,
0146                           std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
0147   fluid_cell_info_ptr = make_unique<FluidCellInfo>();
0148   if (hydro_status != FINISHED) {
0149     throw std::runtime_error("Hydro evolution is not finished ");
0150   }
0151
0152   if (!bulk_info.tau_eta_is_tz) {
0153     Jetscape::real tau = std::sqrt(t * t - z * z);
0154     Jetscape::real eta = 0.5 * (std::log(t + z) - std::log(t - z));
0155     *fluid_cell_info_ptr = bulk_info.get(tau, x, y, eta);
0156   } else {
0157     *fluid_cell_info_ptr = bulk_info.get(t, x, y, z);
0158   }
0159 }