Skip to content

Commit

Permalink
fixed the source Tmunu class
Browse files Browse the repository at this point in the history
  • Loading branch information
@chunshen1987
chunshen1987 committed Sep 23, 2024
1 parent 90f0ede commit 9e5810c
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 39 deletions.
74 changes: 43 additions & 31 deletions src/hydro_source_Tmunu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,26 +7,28 @@
#include <cstdlib>
#include <vector>
#include <array>


using Util::hbarc;

// 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::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::cerr << "Error: Cannot open the IP-Glasma file: "
<< IPG_filename << std::endl;
std::exit(1);
}

Expand All @@ -42,14 +44,14 @@ void HydroSourceTmunu::readIPGevent() {
>> 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::cerr << "Error: Invalid grid size nx = "
<< nx << ", ny = " << ny << std::endl;
std::exit(1);
}

Expand All @@ -70,28 +72,31 @@ void HydroSourceTmunu::readIPGevent() {
double pixx, pixy, pixeta, piyy, piyeta, pietaeta;

if (profile.eof()) {
std::cerr << "Error: Unexpected end of file while reading data." << std::endl;
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
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;
>> 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; // Add this line
DATA.eta_size = DATA.delta_eta * DATA.neta;
}


// Adjust ueta and shear tensor components
ueta *= tau0;
pitaueta *= tau0 * DATA.sFactor;
Expand All @@ -100,35 +105,44 @@ void HydroSourceTmunu::readIPGevent() {
piyeta *= tau0 * DATA.sFactor;

// Store the values in the containers
energy_density_[ix][iy] = e * DATA.sFactor / hbarc;
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};
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 << "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;
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;
}

Expand All @@ -149,16 +163,14 @@ void HydroSourceTmunu::get_hydro_energy_source(
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;

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;

// Unit conversion factor
const double prefactors = 1.0 / (DATA.delta_tau * Util::hbarc);
const double prefactors = 1.0 / dtau;
for (int i = 0; i < 4; ++i) {
j_mu[i] *= prefactors;
j_mu[i] *= prefactors; // 1/fm^5
}
}

3 changes: 2 additions & 1 deletion src/hydro_source_Tmunu.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,13 @@ class HydroSourceTmunu : public HydroSourceBase { // Inherit from HydroSourceBa
// 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;
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_;
Expand Down
13 changes: 6 additions & 7 deletions src/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,7 @@ void Init::InitArena(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr,
music_message << ", dx=" << DATA.delta_x << " fm, dy="
<< DATA.delta_y << " fm.";
music_message.flush("info");
}


else if (DATA.Initial_profile == 1) {
} else if (DATA.Initial_profile == 1) {
music_message << "Using Initial_profile=" << DATA.Initial_profile;
DATA.nx = 2;
DATA.ny = 2;
Expand Down Expand Up @@ -120,7 +117,8 @@ void Init::InitArena(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr,
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");
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()
- DATA.delta_tau);
Expand Down Expand Up @@ -162,7 +160,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");
}
}

// initialize arena
arenaFieldsPrev.resizeFields(DATA.nx, DATA.ny, DATA.neta);
Expand Down Expand Up @@ -251,7 +249,8 @@ void Init::InitTJb(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr) {
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");
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);
Expand Down

0 comments on commit 9e5810c

Please sign in to comment.