Merge branch 'anar_dev' into mergeTemp

.gitignore

@@ -17,3 +17,5 @@ MUSIChydro
 *.dylib
 *.so
 music_input*
+.vscode/
+input.default

src/CMakeLists.txt

@@ -10,6 +10,7 @@ set (SOURCES
     hydro_source_strings.cpp
     hydro_source_ampt.cpp
     hydro_source_TATB.cpp
+    hydro_source_Tmunu.cpp 
     pretty_ostream.cpp
     grid_info.cpp
     util.cpp

src/hydro_source_Tmunu.cpp

@@ -0,0 +1,179 @@
+#include "hydro_source_Tmunu.h"
+
+#include <array>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <vector>
+
+#include "util.h"
+
+// Constructor
+HydroSourceTmunu::HydroSourceTmunu(InitData &DATA_in) : DATA(DATA_in) {
+    // Read the IP-Glasma event file
+    readIPGevent();
+    tau_source_ = DATA.tau0;
+    DATA.tau0 = DATA.tau0 - DATA.delta_tau;
+    set_source_tau_min(tau_source_);
+    set_source_tau_max(tau_source_);
+}
+
+// Function to read the IP-Glasma event file and initialize data
+void HydroSourceTmunu::readIPGevent() {
+    // Open file
+    std::string IPG_filename = DATA.initName;
+    std::ifstream profile(IPG_filename);
+
+    if (!profile.is_open()) {
+        std::cerr << "Error: Cannot open the IP-Glasma file: " << IPG_filename
+                  << std::endl;
+        std::exit(1);
+    }
+
+    // Read the information line
+    std::string dummy;
+    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;
+
+    // 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::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;
+                DATA.eta_size = DATA.delta_eta * DATA.neta;
+                music_message << "eta_size=" << DATA.eta_size
+                              << ", x_size=" << DATA.x_size
+                              << ", y_size=" << DATA.y_size;
+                music_message.flush("info");
+            }
+
+            // 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 / Util::hbarc;  // 1/fm^4
+            velocity_[ix][iy] = {utau, ux, uy, ueta};
+            shear_tensor_[ix][iy] = {pitautau, pitaux,  pitauy, pitaueta,
+                                     pixx,     pixy,    pixeta, piyy,
+                                     piyeta,   pietaeta};  // 1/fm^4
+        }
+    }
+    profile.close();
+    std::cout << "neta = " << DATA.neta << ", eta_size = " << DATA.eta_size
+              << ", delta_eta = " << DATA.delta_eta << std::endl;
+    std::cout << "nx = " << DATA.nx << ", x_size = " << DATA.x_size
+              << ", delta_x = " << DATA.delta_x << std::endl;
+    std::cout << "ny = " << DATA.ny << ", x_size = " << DATA.y_size
+              << ", delta_y = " << DATA.delta_y << std::endl;
+}
+
+// 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 {
+    const double dtau = DATA.delta_tau;
+    j_mu = {0};
+
+    if (std::abs((tau - tau_source_)) > 0.5 * dtau) return;
+
+    // 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;
+    }
+
+    // Access data
+    double epsilon = energy_density_[ix][iy];
+    const auto &Wmunu = shear_tensor_[ix][iy];
+    const auto &u = velocity_[ix][iy];
+
+    // Extract individual elements
+    double u0 = u[0];
+    double u1 = u[1];
+    double u2 = u[2];
+    double u3 = u[3];
+
+    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) + Wmunu0;
+    j_mu[1] = epsilon * u0 * u1 + (epsilon / 3.0) * u0 * u1 + Wmunu1;
+    j_mu[2] = epsilon * u0 * u2 + (epsilon / 3.0) * u0 * u2 + Wmunu2;
+    j_mu[3] = epsilon * u0 * u3 + (epsilon / 3.0) * u0 * u3 + Wmunu3;
+
+    const double prefactors = 1.0 / dtau;
+    for (int i = 0; i < 4; ++i) {
+        j_mu[i] *= prefactors;  // 1/fm^5
+    }
+}

src/hydro_source_Tmunu.h

@@ -0,0 +1,37 @@
+// hydro_source_Tmunu.h
+
+#ifndef HYDRO_SOURCE_TMUNU_H
+#define HYDRO_SOURCE_TMUNU_H
+
+#include <array>
+#include <vector>
+
+#include "data_struct.h"
+#include "hydro_source_base.h"  // Include the base class header
+
+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;
+
+  private:
+    // Add your private members and functions here
+    InitData &DATA;
+
+    double tau_source_;
+    // 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();
+};
+
+#endif  // HYDRO_SOURCE_TMUNU_H

src/init.cpp

@@ -123,6 +123,9 @@ void Init::InitArena(
         DATA.tau0 = std::max(DATA.tau0, tau_overlap) - DATA.delta_tau;
         music_message << "tau0 = " << DATA.tau0 << " fm/c.";
         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) {
         DATA.tau0 =
             (hydro_source_terms_ptr.lock()->get_source_tau_min()
@@ -244,6 +247,13 @@ void Init::InitTJb(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr) {
         initial_MCGlb_with_rhob(arenaFieldsPrev, arenaFieldsCurr);
     } else if (DATA.Initial_profile == 112 || DATA.Initial_profile == 113) {
         music_message.info("Initialize hydro with source terms from TA and TB");
+#pragma omp parallel for
+        for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
+            initial_with_zero_XY(ieta, arenaFieldsPrev, arenaFieldsCurr);
+        }
+    } else if (DATA.Initial_profile == 17) {
+        music_message.info(
+            "Initialize hydro with source terms from pre-equilibrium model");
 #pragma omp parallel for
         for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
             initial_with_zero_XY(ieta, arenaFieldsPrev, arenaFieldsCurr);

src/music.cpp

@@ -13,6 +13,7 @@
 #include "dissipative.h"
 #include "evolve.h"
 #include "hydro_source_TATB.h"
+#include "hydro_source_Tmunu.h"
 #include "hydro_source_ampt.h"
 #include "hydro_source_strings.h"
 #include "init.h"
@@ -66,6 +67,10 @@ void MUSIC::generate_hydro_source_terms() {
         auto hydro_source_ptr =
             std::shared_ptr<HydroSourceTATB>(new HydroSourceTATB(DATA));
         add_hydro_source_terms(hydro_source_ptr);
+    } else if (DATA.Initial_profile == 17) {
+        auto hydro_source_ptr =
+            std::shared_ptr<HydroSourceTmunu>(new HydroSourceTmunu(DATA));
+        add_hydro_source_terms(hydro_source_ptr);
     }
 }