From 90f0edeead45c87ebec7d3a84c26fb16a2731e12 Mon Sep 17 00:00:00 2001 From: Anar Akbarov Date: Mon, 23 Sep 2024 11:59:07 -0400 Subject: [PATCH] Fix the reading of the file --- src/hydro_source_Tmunu.cpp | 236 +++++++++++++++++-------------------- src/hydro_source_Tmunu.h | 47 ++++---- src/init.cpp | 11 +- 3 files changed, 134 insertions(+), 160 deletions(-) diff --git a/src/hydro_source_Tmunu.cpp b/src/hydro_source_Tmunu.cpp index 08122e89..0cc26e60 100644 --- a/src/hydro_source_Tmunu.cpp +++ b/src/hydro_source_Tmunu.cpp @@ -7,7 +7,7 @@ #include #include #include -#include // Include for OpenMP + using Util::hbarc; @@ -22,167 +22,143 @@ HydroSourceTmunu::HydroSourceTmunu(InitData &DATA_in) void HydroSourceTmunu::readIPGevent() { // Open file std::string IPG_filename = DATA.initName; - std::ifstream IPGinputFile(IPG_filename); + std::ifstream profile(IPG_filename); - if (!IPGinputFile.is_open()) { + + if (!profile.is_open()) { std::cerr << "Error: Cannot open the IP-Glasma file: " << IPG_filename << std::endl; std::exit(1); } - // Step1: Initialize containers for energy density, velocity, and shear tensor - energy_density_ = std::vector>>( - DATA.nx, std::vector>(DATA.ny, std::vector(DATA.neta, 0.0))); - - velocity_ = std::vector>>>( - DATA.nx, std::vector>>( - DATA.ny, std::vector>(DATA.neta))); - - shear_tensor_ = std::vector>>>( - DATA.nx, std::vector>>( - DATA.ny, std::vector>(DATA.neta))); - - std::string line; - double tau0 = DATA.tau0; - std::ifstream profile(DATA.initName.c_str()); // Fixed the ifstream - + // Read the information line std::string dummy; - // read the information line std::getline(profile, dummy); + std::istringstream ss1(dummy); + + int nx, ny, neta; + double deta, dx, dy, dummy2; + // read the first line with general info + ss1 >> dummy >> dummy >> dummy2 + >> dummy >> neta >> dummy >> nx >> dummy >> ny + >> dummy >> deta >> dummy >> dx >> dummy >> dy; + + + DATA.nx = nx; + DATA.ny = ny; + DATA.delta_x = dx; + DATA.delta_y = dy; + + if (nx <= 0 || ny <= 0) { + std::cerr << "Error: Invalid grid size nx = " << nx << ", ny = " << ny << std::endl; + std::exit(1); + } + + // Initialize containers + energy_density_.resize(nx, std::vector(ny, 0.0)); + velocity_.resize(nx, std::vector>(ny)); + shear_tensor_.resize(nx, std::vector>(ny)); + + double tau0 = DATA.tau0; - const int nx = DATA.nx; // Changed to correct variable from DATA - const int ny = DATA.ny; // Changed to correct variable from DATA - std::vector temp_profile_ed(nx * ny, 0.0); - std::vector temp_profile_utau(nx * ny, 0.0); - std::vector temp_profile_ux(nx * ny, 0.0); - std::vector temp_profile_uy(nx * ny, 0.0); - std::vector temp_profile_ueta(nx * ny, 0.0); - std::vector temp_profile_pitautau(nx * ny, 0.0); - std::vector temp_profile_pitaux(nx * ny, 0.0); - std::vector temp_profile_pitauy(nx * ny, 0.0); - std::vector temp_profile_pitaueta(nx * ny, 0.0); - std::vector temp_profile_pixx(nx * ny, 0.0); - std::vector temp_profile_pixy(nx * ny, 0.0); - std::vector temp_profile_pixeta(nx * ny, 0.0); - std::vector temp_profile_piyy(nx * ny, 0.0); - std::vector temp_profile_piyeta(nx * ny, 0.0); - std::vector temp_profile_pietaeta(nx * ny, 0.0); - - // Step2: Reading IP-Glasma data - int ix, iy, ieta; - - double e, utau, ux, uy, ueta; - double density = 0, dummy2 = 0, dummy3 = 0; - double pixx, pixy, pixeta, piyy, piyeta, pietaeta, pitautau, pitaux, pitauy, pitaueta; - - for (ix = 0; ix < nx; ix++) { - for (iy = 0; iy < ny; iy++) { - int idx = iy + ix * ny; + // Parallelize the loops if desired + //#pragma omp parallel for collapse(2) + for (int ix = 0; ix < nx; ix++) { + for (int iy = 0; iy < ny; iy++) { + double dummy1, dummy2, dummy3; + double e, utau, ux, uy, ueta; + double pitautau, pitaux, pitauy, pitaueta; + double pixx, pixy, pixeta, piyy, piyeta, pietaeta; + + if (profile.eof()) { + std::cerr << "Error: Unexpected end of file while reading data." << std::endl; + std::exit(1); + } + std::getline(profile, dummy); - std::stringstream ss(dummy); // Fixed the issue with ss and iss - ss >> ix >> iy >> ieta >> e >> utau >> ux >> uy >> ueta - >> pitautau >> pitaux >> pitauy >> pitaueta - >> pixx >> pixy >> pixeta >> piyy >> piyeta >> pietaeta; - - ueta = ueta * tau0; - temp_profile_ed[idx] = density * DATA.sFactor / hbarc; // 1/fm^4 - temp_profile_ux[idx] = ux; - temp_profile_uy[idx] = uy; - temp_profile_ueta[idx] = ueta; - temp_profile_utau[idx] = sqrt(1. + ux * ux + uy * uy + ueta * ueta); - - // Viscous components already in 1/fm^4 from the IP-Glasma output - double visFactor = DATA.sFactor * DATA.preEqVisFactor; - temp_profile_pixx[idx] = pixx * visFactor; - temp_profile_pixy[idx] = pixy * visFactor; - temp_profile_pixeta[idx] = pixeta * tau0 * visFactor; - temp_profile_piyy[idx] = piyy * visFactor; - temp_profile_piyeta[idx] = piyeta * tau0 * visFactor; - - utau = temp_profile_utau[idx]; - temp_profile_pietaeta[idx] = ( - (2. * (ux * uy * temp_profile_pixy[idx] - + ux * ueta * temp_profile_pixeta[idx] - + uy * ueta * temp_profile_piyeta[idx]) - - (utau * utau - ux * ux) * temp_profile_pixx[idx] - - (utau * utau - uy * uy) * temp_profile_piyy[idx]) - / (utau * utau - ueta * ueta)); - temp_profile_pitaux[idx] = (1. / utau - * (temp_profile_pixx[idx] * ux - + temp_profile_pixy[idx] * uy - + temp_profile_pixeta[idx] * ueta)); - temp_profile_pitauy[idx] = (1. / utau - * (temp_profile_pixy[idx] * ux - + temp_profile_piyy[idx] * uy - + temp_profile_piyeta[idx] * ueta)); - temp_profile_pitaueta[idx] = (1. / utau - * (temp_profile_pixeta[idx] * ux - + temp_profile_piyeta[idx] * uy - + temp_profile_pietaeta[idx] * ueta)); - temp_profile_pitautau[idx] = (1. / utau - * (temp_profile_pitaux[idx] * ux - + temp_profile_pitauy[idx] * uy - + temp_profile_pitaueta[idx] * ueta)); + std::istringstream ss(dummy); + // Read data from the file + if (!(ss >> dummy1 >> dummy2 >> dummy3 >> e >> utau >> ux >> uy >> ueta + >> pitautau >> pitaux >> pitauy >> pitaueta + >> pixx >> pixy >> pixeta >> piyy >> piyeta >> pietaeta)) { + std::cerr << "Error: Failed to parse data line: " << dummy << std::endl; + std::exit(1); + } + // Set x_size and y_size based on the data if (ix == 0 && iy == 0) { DATA.x_size = -dummy2 * 2; DATA.y_size = -dummy3 * 2; - if (omp_get_thread_num() == 0) { - music_message << "eta_size=" << DATA.eta_size - << ", x_size=" << DATA.x_size - << ", y_size=" << DATA.y_size; - music_message.flush("info"); - } + DATA.eta_size = DATA.delta_eta * DATA.neta; // Add this line } - } - // Store the data in the containers - double tau = DATA.tau0; // Fixed tau declaration - shear_tensor_[ix][iy][ieta] = {pitautau, pitaux, pitauy, pitaueta * tau * DATA.sFactor, - pixx * DATA.sFactor, pixy * DATA.sFactor, pixeta * DATA.sFactor, piyy * DATA.sFactor, piyeta * DATA.sFactor, pietaeta}; // multiply by DATA.sFactor + + // Adjust ueta and shear tensor components + ueta *= tau0; + pitaueta *= tau0 * DATA.sFactor; + pietaeta *= tau0 * tau0 * DATA.sFactor; + pixeta *= tau0 * DATA.sFactor; + piyeta *= tau0 * DATA.sFactor; + + // Store the values in the containers + energy_density_[ix][iy] = e * DATA.sFactor / hbarc; + velocity_[ix][iy] = {utau, ux, uy, ueta}; + shear_tensor_[ix][iy] = {pitautau, pitaux, pitauy, pitaueta, + pixx, pixy, pixeta, piyy, piyeta, pietaeta}; + } } - IPGinputFile.close(); + profile.close(); } -// Function to calculate energy-momentum source term +// Function to calculate energy-momentum source terms void HydroSourceTmunu::get_hydro_energy_source( const double tau, const double x, const double y, const double eta_s, const std::array& u_mu, std::array& j_mu) const { j_mu = {0}; - - // Grid indices from spatial coordinates - int ix = static_cast((x + DATA.x_size / 2.0) / DATA.delta_x); - int iy = static_cast((y + DATA.y_size / 2.0) / DATA.delta_y); - int ieta = static_cast((eta_s + DATA.eta_size / 2.0) / DATA.delta_eta); - - // Bounds for grid indices - if (ix < 0 || ix >= DATA.nx || iy < 0 || iy >= DATA.ny || ieta < 0 || ieta >= DATA.neta) { - std::cerr << "Out of bounds grid access in HydroSourceTmunu::get_hydro_energy_source" << std::endl; - std::exit(1); + + // Compute indices from positions + const int ix = static_cast((x + DATA.x_size/2.)/DATA.delta_x + 0.1); + const int iy = static_cast((y + DATA.y_size/2.)/DATA.delta_y + 0.1); + + + // Bounds checking + if (ix < 0 || ix >= DATA.nx || iy < 0 || iy >= DATA.ny) { + std::cerr << "Error: Index out of bounds in get_hydro_energy_source: ix = " << ix << ", iy = " << iy << std::endl; + std::cerr << "x = " << x << ", y = " << y << std::endl; + std::cerr << "x_size = " << DATA.x_size << ", y_size = " << DATA.y_size << std::endl; + std::cerr << "delta_x = " << DATA.delta_x << ", delta_y = " << DATA.delta_y << std::endl; + return; } - // Read the energy density, velocity, and viscous tensor from the IP-Glasma data - double epsilon = energy_density_[ix][iy][ieta]; // Energy density e - const auto &Wmunu = shear_tensor_[ix][iy][ieta]; // Shear viscous tensor pi^{mu nu} - const auto &u = velocity_[ix][iy][ieta]; // Fluid velocity u^mu + // Access data + double epsilon = energy_density_[ix][iy]; + const auto &Wmunu = shear_tensor_[ix][iy]; + const auto &u = velocity_[ix][iy]; - // Compute j^0 - j_mu[0] = epsilon * (u[0] * u[0]) + (epsilon / 3) * (1 + u[0] * u[0]) + Wmunu[0]; // \pi^{\tau \tau} + // Extract individual elements + double u0 = u[0]; + double u1 = u[1]; + double u2 = u[2]; + double u3 = u[3]; - // Compute spatial components - j_mu[1] = epsilon * u[0] * u[1] + (epsilon / 3) * u[0] * u[1] + Wmunu[1]; // \pi^{\tau x} - j_mu[2] = epsilon * u[0] * u[2] + (epsilon / 3) * u[0] * u[2] + Wmunu[2]; // \pi^{\tau y} - j_mu[3] = epsilon * u[0] * u[3] + (epsilon / 3) * u[0] * u[3] + Wmunu[3]; // \pi^{\tau \eta} + double Wmunu0 = Wmunu[0]; + double Wmunu1 = Wmunu[1]; + double Wmunu2 = Wmunu[2]; + double Wmunu3 = Wmunu[3]; + // Compute j^0 and j^i (energy-momentum source term) + j_mu[0] = epsilon * u0 * u0 - (epsilon / 3.0) * (1.0 - u0 * u0) + 0*Wmunu0; + j_mu[1] = epsilon * u0 * u1 + (epsilon / 3.0) * u0 * u1 + 0*Wmunu1; + j_mu[2] = epsilon * u0 * u2 + (epsilon / 3.0) * u0 * u2 + 0*Wmunu2; + j_mu[3] = epsilon * u0 * u3 + (epsilon / 3.0) * u0 * u3 + 0*Wmunu3; - // Unit conversion factor const double prefactors = 1.0 / (DATA.delta_tau * Util::hbarc); - j_mu[0] *= prefactors; - j_mu[1] *= prefactors; - j_mu[2] *= prefactors; - j_mu[3] *= prefactors; + for (int i = 0; i < 4; ++i) { + j_mu[i] *= prefactors; + } } + diff --git a/src/hydro_source_Tmunu.h b/src/hydro_source_Tmunu.h index c4c967e9..685f16cb 100644 --- a/src/hydro_source_Tmunu.h +++ b/src/hydro_source_Tmunu.h @@ -1,31 +1,34 @@ -#ifndef SRC_HYDRO_SOURCE_TMUNU_H_ -#define SRC_HYDRO_SOURCE_TMUNU_H_ +// hydro_source_Tmunu.h -#include "data.h" -#include "hydro_source_base.h" -#include "eos.h" -#include "cell.h" +#ifndef HYDRO_SOURCE_TMUNU_H +#define HYDRO_SOURCE_TMUNU_H + +#include "hydro_source_base.h" // Include the base class header +#include "data_struct.h" #include #include -#include -#include -#include // for exit function -class HydroSourceTmunu : public HydroSourceBase { - private: +class HydroSourceTmunu : public HydroSourceBase { // Inherit from HydroSourceBase +public: + // Constructor + HydroSourceTmunu(InitData &DATA_in); + + // Override the virtual function from HydroSourceBase + void get_hydro_energy_source( + const double tau, const double x, const double y, const double eta_s, + const std::array& u_mu, std::array& j_mu) const override; + +private: + // Add your private members and functions here InitData &DATA; - std::vector>> energy_density_; - std::vector>>> velocity_; - std::vector>>> shear_tensor_; + // Data containers + std::vector> energy_density_; + std::vector>> velocity_; + std::vector>> shear_tensor_; + // Function to read the IP-Glasma event file and initialize data void readIPGevent(); - - public: - HydroSourceTmunu() = delete; - HydroSourceTmunu(InitData &DATA_in); - //~HydroSourceTmunu() = default; - void get_hydro_energy_source(const double tau, const double x, const double y, const double eta_s, - const std::array& u_mu, std::array& j_mu) const; }; -#endif + +#endif // HYDRO_SOURCE_TMUNU_H diff --git a/src/init.cpp b/src/init.cpp index c1058b89..c6a55dd5 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -80,7 +80,7 @@ void Init::InitArena(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr, << " fm, dy=" << DATA.delta_y << " fm."; music_message.flush("info"); } else if ( DATA.Initial_profile == 9 || DATA.Initial_profile == 91 - || DATA.Initial_profile == 92 || DATA.Initial_profile == 93) { + || DATA.Initial_profile == 92 || DATA.Initial_profile == 93 ) { music_message.info(DATA.initName); ifstream profile(DATA.initName.c_str()); if (!profile.is_open()) { @@ -121,9 +121,7 @@ void Init::InitArena(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr, music_message.flush("info"); } else if (DATA.Initial_profile == 17) { music_message.info("Initialize hydro with source terms from pre-equilibrium model"); - } - - else if (DATA.Initial_profile == 13 || DATA.Initial_profile == 131) { + } else if (DATA.Initial_profile == 13 || DATA.Initial_profile == 131) { DATA.tau0 = (hydro_source_terms_ptr.lock()->get_source_tau_min() - DATA.delta_tau); DATA.tau0 = static_cast(DATA.tau0/0.02)*0.02; @@ -164,10 +162,7 @@ void Init::InitArena(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr, music_message << "dx = " << DATA.delta_x << " fm, dy = " << DATA.delta_y << " fm, deta = " << DATA.delta_eta; music_message.flush("info"); - } else if (DATA.Initial_profile == 17){ - music_message.info ("Initialize hydro with source terms from pre-equilibrium model"); - - } + } // initialize arena arenaFieldsPrev.resizeFields(DATA.nx, DATA.ny, DATA.neta);