Fix the reading of the file

src/hydro_source_Tmunu.cpp

@@ -7,7 +7,7 @@
 #include <cstdlib>
 #include <vector>
 #include <array>
-#include <omp.h>  // 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<std::vector<std::vector<double>>>(
-        DATA.nx, std::vector<std::vector<double>>(DATA.ny, std::vector<double>(DATA.neta, 0.0)));
-
-    velocity_ = std::vector<std::vector<std::vector<std::array<double, 4>>>>(
-        DATA.nx, std::vector<std::vector<std::array<double, 4>>>(
-            DATA.ny, std::vector<std::array<double, 4>>(DATA.neta)));
-
-    shear_tensor_ = std::vector<std::vector<std::vector<std::array<double, 10>>>>(
-        DATA.nx, std::vector<std::vector<std::array<double, 10>>>(
-            DATA.ny, std::vector<std::array<double, 10>>(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<double>(ny, 0.0));
+    velocity_.resize(nx, std::vector<std::array<double, 4>>(ny));
+    shear_tensor_.resize(nx, std::vector<std::array<double, 10>>(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<double> temp_profile_ed(nx * ny, 0.0);
-    std::vector<double> temp_profile_utau(nx * ny, 0.0);
-    std::vector<double> temp_profile_ux(nx * ny, 0.0);
-    std::vector<double> temp_profile_uy(nx * ny, 0.0);
-    std::vector<double> temp_profile_ueta(nx * ny, 0.0);
-    std::vector<double> temp_profile_pitautau(nx * ny, 0.0);
-    std::vector<double> temp_profile_pitaux(nx * ny, 0.0);
-    std::vector<double> temp_profile_pitauy(nx * ny, 0.0);
-    std::vector<double> temp_profile_pitaueta(nx * ny, 0.0);
-    std::vector<double> temp_profile_pixx(nx * ny, 0.0);
-    std::vector<double> temp_profile_pixy(nx * ny, 0.0);
-    std::vector<double> temp_profile_pixeta(nx * ny, 0.0);
-    std::vector<double> temp_profile_piyy(nx * ny, 0.0);
-    std::vector<double> temp_profile_piyeta(nx * ny, 0.0);
-    std::vector<double> 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<double, 4>& u_mu, std::array<double, 4>& j_mu) const {
 
     j_mu = {0};
 
-
-    // Grid indices from spatial coordinates
-    int ix = static_cast<int>((x + DATA.x_size / 2.0) / DATA.delta_x);
-    int iy = static_cast<int>((y + DATA.y_size / 2.0) / DATA.delta_y);
-    int ieta = static_cast<int>((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<int>((x + DATA.x_size/2.)/DATA.delta_x + 0.1);
+    const int iy = static_cast<int>((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;
+    }
 }
+

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 <vector>
 #include <array>
-#include <memory>
-#include <iostream>
-#include <cstdlib> // 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<double, 4>& u_mu, std::array<double, 4>& j_mu) const override;
+
+private:
+    // Add your private members and functions here
     InitData &DATA;
 
-    std::vector<std::vector<std::vector<double>>> energy_density_;
-    std::vector<std::vector<std::vector<std::array<double, 4>>>> velocity_;
-    std::vector<std::vector<std::vector<std::array<double, 10>>>> shear_tensor_;
+    // Data containers
+    std::vector<std::vector<double>> energy_density_;
+    std::vector<std::vector<std::array<double, 4>>> velocity_;
+    std::vector<std::vector<std::array<double, 10>>> 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<double, 4>& u_mu, std::array<double, 4>& j_mu) const;
 };
-#endif
+
+#endif  // HYDRO_SOURCE_TMUNU_H

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<int>(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);