Skip to content

Commit

Permalink
Merge branch 'greg_dev' of https://github.com/MUSIC-fluid/MUSIC into …
Browse files Browse the repository at this point in the history 
…greg_dev
  • Loading branch information
@chunshen1987
chunshen1987 committed Jun 2, 2024
2 parents 4ffa7ca + cd543ee commit c7587a1
Show file tree
Hide file tree
Showing 16 changed files with 97 additions and 108 deletions.
2 changes: 1 addition & 1 deletion src/dissipative.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ using Util::map_1d_idx_to_2d;

Diss::Diss(const EOS &eosIn, const InitData &Data_in) :
DATA(Data_in), eos(eosIn), minmod(Data_in),
transport_coeffs_(eosIn, Data_in) {}
transport_coeffs_(Data_in) {}

/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
/* Dissipative parts */
Expand Down
2 changes: 1 addition & 1 deletion src/eos.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ class EOS {
std::unique_ptr<EOS_base> eos_ptr;

public:
EOS() = default;
EOS() = delete;
EOS(const int eos_id_in);

~EOS() {};
Expand Down
44 changes: 24 additions & 20 deletions src/eos_4D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,26 +187,29 @@ void EOS_4D::FourDLInterp(const std::vector<float> &data,
std::array<float, 5> &ResArr,
bool compute_derivatives) const {
// Constrain the input tilde variables values to the table boundaries.
float T = std::max(Ttilde0, std::min(T_tilde_max, TildeVar[0]));
float mub = std::max(mubtilde0, std::min(mub_tilde_max, TildeVar[1]));
float muq = std::max(muqtilde0, std::min(muq_tilde_max, TildeVar[2]));
float mus = std::max(mustilde0, std::min(mus_tilde_max, TildeVar[3]));

// Calculate the weights associated to the sixteen surrounding point
float indmuT = (T - Ttilde0)/dTtilde;
float indmub = (mub - mubtilde0)/dmubtilde;
float indmuq = (muq - muqtilde0)/dmuqtilde;
float indmus = (mus - mustilde0)/dmustilde;

int iT = static_cast<int>(indmuT); int iT1 = iT + 1;
int ib = static_cast<int>(indmub); int ib1 = ib + 1;
int iq = static_cast<int>(indmuq); int iq1 = iq + 1;
int is = static_cast<int>(indmus); int is1 = is + 1;

float dx = indmuT - iT;
float dy = indmub - ib;
float dz = indmuq - iq;
float dt = indmus - is;
double Tfrac = (TildeVar[0] - Ttilde0)/dTtilde;
double mubfrac = (TildeVar[1] - mubtilde0)/dmubtilde;
double muqfrac = (TildeVar[2] - muqtilde0)/dmuqtilde;
double musfrac = (TildeVar[3] - mustilde0)/dmustilde;

// Calculate the weights associated to the sixteen surrounding point
int iT = static_cast<int>(Tfrac);
iT = std::max(0, std::min(N_T - 2, iT));
int iT1 = iT + 1;
int ib = static_cast<int>(mubfrac);
ib = std::max(0, std::min(N_mub - 2, ib));
int ib1 = ib + 1;
int iq = static_cast<int>(muqfrac);
iq = std::max(0, std::min(N_muq - 2, iq));
int iq1 = iq + 1;
int is = static_cast<int>(musfrac);
is = std::max(0, std::min(N_mus - 2, is));
int is1 = is + 1;

float dx = std::max(0., std::min(1., Tfrac - iT));
float dy = std::max(0., std::min(1., mubfrac - ib));
float dz = std::max(0., std::min(1., muqfrac - iq));
float dt = std::max(0., std::min(1., musfrac - is));

float w0000 = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dt);
float w1111 = dx * dy * dz * dt;
Expand Down Expand Up @@ -290,6 +293,7 @@ void EOS_4D::FourDLInterp(const std::vector<float> &data,
//float dXdTtilde = (tempT2 - tempT1)/dTtilde;
float T1 = Ttilde0 + iT*dTtilde;
float T2 = T1 + dTtilde;
float T = std::max(Ttilde0, std::min(T_tilde_max, TildeVar[0]));
float dXdTtilde = (
5*(tempT2 - tempT1)/(pow(T2, 5) - pow(T1, 5))*pow(T, 4));
// to test for speed
Expand Down
1 change: 0 additions & 1 deletion src/eos_4D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@ class EOS_4D : public EOS_base {
//double Ttilde, mubtilde, muqtilde, mustilde;

// useful constants
const int Nf = 3;
const double alphaNf = (8/45.0 + 7/60.0*3.0)*M_PI*M_PI;
const double OneoveralphaNf = 1./alphaNf;

Expand Down
26 changes: 8 additions & 18 deletions src/freeze_pseudo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ void Freeze::ComputeParticleSpectrum_pseudo_improved(InitData *DATA,
int iphimax = DATA->phi_steps; // number of points
// (phi=2pi equal to phi=0)
double deltaphi = 2*M_PI/iphimax;

// Reuse rapidity variables (Need to reuse variable y
// so resonance decay routine can be used as is.
// Might as well misuse ymax and deltaY too)
Expand All @@ -163,11 +163,11 @@ void Freeze::ComputeParticleSpectrum_pseudo_improved(InitData *DATA,
<< particleList[j].name << "("
<< particleList[j].number << ") ... ";
music_message.flush("info");

particleList[j].ny = DATA->pseudo_steps + 1;
particleList[j].npt = iptmax;
particleList[j].nphi = iphimax;

// set particle properties
double m = particleList[j].mass;
int deg = particleList[j].degeneracy;
Expand All @@ -183,23 +183,13 @@ void Freeze::ComputeParticleSpectrum_pseudo_improved(InitData *DATA,
}

// open files to write
FILE *d_file;
char *buf = new char[10];
sprintf (buf, "%d", 0);

char *specString=new char[30];
strcpy(specString, "yptphiSpectra");
strcat(specString, buf);
strcat(specString, ".dat");
char* s_name = specString;
std::string s_name = "yptphiSpectra0.dat";
FILE *s_file;
s_file = fopen(s_name, "w");

delete[] specString;
delete[] buf;
s_file = fopen(s_name.c_str(), "w");

// -----------------------------------------------------------------------
// write information in the particle information file
FILE *d_file;
const char* d_name = "particleInformation.dat";
d_file = fopen(d_name, "a");
fprintf(d_file, "%d %e %d %e %e %d %d \n",
Expand All @@ -222,7 +212,7 @@ void Freeze::ComputeParticleSpectrum_pseudo_improved(InitData *DATA,
pt_array[ipt] = pt;
particleList[j].pt[ipt] = pt;
}

double *bulk_deltaf_coeffs = new double[3];
if (DATA_ptr->turn_on_bulk == 1 && DATA_ptr->include_deltaf_bulk == 1) {
if (bulk_deltaf_kind == 0) {
Expand All @@ -231,7 +221,7 @@ void Freeze::ComputeParticleSpectrum_pseudo_improved(InitData *DATA,
bulk_deltaf_coeffs = new double[2];
}
}

double alpha = 0.0;
// main loop begins ...
// store E dN/d^3p as function of phi,
Expand Down
23 changes: 7 additions & 16 deletions src/grid_info.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1677,8 +1677,6 @@ void Cell_info::output_vorticity_distribution(
VorticityVec omega_T = {0.0};
double T_avg = 0.0;
double muB_avg = 0.0;
double muS_avg = 0.0;
double muQ_avg = 0.0;
double weight = 0.0;
for (int ieta = 0; ieta < arena_curr.nEta(); ieta++) {
int fieldIdx = arena_curr.getFieldIdx(ix, iy, ieta);
Expand All @@ -1689,23 +1687,16 @@ void Cell_info::output_vorticity_distribution(
const double e_local = arena_curr.e_[fieldIdx];
if (e_local < 0.1) continue;

const double rhob_local = arena_curr.rhob_[fieldIdx];
const double rhoq_local = arena_curr.rhoq_[fieldIdx];
const double rhos_local = arena_curr.rhos_[fieldIdx];

const double T_local = eos.get_temperature(e_local, rhob_local,
rhoq_local, rhos_local);
const double muB_local = eos.get_muB(e_local, rhob_local,
rhoq_local, rhos_local);
const double muQ_local = eos.get_muQ(e_local, rhob_local,
rhoq_local, rhos_local);
const double muS_local = eos.get_muS(e_local, rhob_local,
rhoq_local, rhos_local);
const double rhob_local = arena_curr.rhob_[fieldIdx];
const double rhoq_local = arena_curr.rhoq_[fieldIdx];
const double rhos_local = arena_curr.rhos_[fieldIdx];
const double T_local = eos.get_temperature(
e_local, rhob_local, rhoq_local, rhos_local);
const double muB_local = eos.get_muB(
e_local, rhob_local, rhoq_local, rhos_local);

T_avg += e_local*T_local*hbarc;
muB_avg += e_local*muB_local*hbarc;
muQ_avg += e_local*muQ_local*hbarc;
muS_avg += e_local*muS_local*hbarc;

VorticityVec omega_local_1, omega_local_2;
VorticityVec omega_local_3, omega_local_4;
Expand Down
2 changes: 1 addition & 1 deletion src/hydro_source_TATB.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ class HydroSourceTATB : public HydroSourceBase {
std::vector<std::vector<double>> profile_TB;

public:
HydroSourceTATB() = default;
HydroSourceTATB() = delete;
HydroSourceTATB(InitData &DATA_in);
~HydroSourceTATB();

Expand Down
6 changes: 3 additions & 3 deletions src/hydro_source_ampt.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,13 @@ class HydroSourceAMPT : public HydroSourceBase {
std::vector<std::shared_ptr<parton>> parton_list_current_tau;

public:
HydroSourceAMPT() = default;
HydroSourceAMPT() = delete;
HydroSourceAMPT(const InitData &DATA_in);
~HydroSourceAMPT();

//! This function reads in the partons information from the AMPT model
void read_in_AMPT_partons();

//! this function returns the energy source term J^\mu at a given point
//! (tau, x, y, eta_s)
void get_hydro_energy_source(
Expand Down
2 changes: 1 addition & 1 deletion src/hydro_source_strings.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ class HydroSourceStrings : public HydroSourceBase {
std::vector<std::shared_ptr<QCD_string>> QCD_strings_electric_list_current_tau;

public:
HydroSourceStrings() = default;
HydroSourceStrings() = delete;
HydroSourceStrings(InitData &DATA_in);
~HydroSourceStrings();

Expand Down
47 changes: 27 additions & 20 deletions src/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,21 +197,26 @@ void Init::InitTJb(Fields &arenaFieldsPrev, Fields &arenaFieldsCurr) {
initial_Gubser_XY(ieta, arenaFieldsPrev, arenaFieldsCurr);
}

if (DATA.turn_on_QS == 1){
if(DATA.use_rhoQS_to_rhoB_ratios == 1){
music_message.info("Electric charge and strangeness densities initialized with");
music_message << "rhoQ = " << DATA.ratio_q << " rhoB";
music_message << " rhoS = " << DATA.ratio_s << " rhoB";
music_message.flush("info");
for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
initial_QS_with_rhob_ratios_XY(ieta, arenaFieldsPrev, arenaFieldsCurr);}
}
else{
music_message.info("Initialize electric charge and strangeness densities to zero");
music_message.flush("info");
for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
initial_QS_with_zero_XY(ieta, arenaFieldsPrev, arenaFieldsCurr);}
}
if (DATA.turn_on_QS == 1) {
if (DATA.use_rhoQS_to_rhoB_ratios == 1) {
music_message << "Electric charge and strangeness densities "
<< "initialized with"
<< "rhoQ = " << DATA.ratio_q << " rhoB"
<< " rhoS = " << DATA.ratio_s << " rhoB";
music_message.flush("info");
for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
initial_QS_with_rhob_ratios_XY(
ieta, arenaFieldsPrev, arenaFieldsCurr);
}
} else {
music_message << "Initialize electric charge and strangeness "
<< "densities to zero";
music_message.flush("info");
for (int ieta = 0; ieta < arenaFieldsCurr.nEta(); ieta++) {
initial_QS_with_zero_XY(ieta, arenaFieldsPrev,
arenaFieldsCurr);
}
}
}
} else if (DATA.Initial_profile == 1) {
// code test in 1+1 D vs Monnai's results
Expand Down Expand Up @@ -1255,9 +1260,10 @@ void Init::output_initial_density_profiles(Fields &arena) {
of << "# x(fm) y(fm) eta ed(GeV/fm^3)";
if (DATA.turn_on_rhob == 1){
of << " rhob(1/fm^3)";}
if(DATA.turn_on_QS){
of << " rhoq(1/fm^3)";
of << " rhos(1/fm^3)";}
if (DATA.turn_on_QS) {
of << " rhoq(1/fm^3)";
of << " rhos(1/fm^3)";
}
of << std::endl;
for (int ieta = 0; ieta < arena.nEta(); ieta++) {
double eta_local = (DATA.delta_eta)*ieta - (DATA.eta_size)/2.0;
Expand All @@ -1270,8 +1276,9 @@ void Init::output_initial_density_profiles(Fields &arena) {
<< x_local << " " << y_local << " "
<< eta_local << " " << arena.e_[fieldIdx]*hbarc;
if (DATA.turn_on_rhob == 1) {
of << " " << arena.rhob_[fieldIdx];}
if(DATA.turn_on_QS == 1){
of << " " << arena.rhob_[fieldIdx];
}
if (DATA.turn_on_QS == 1) {
of << " " << arena.rhoq_[fieldIdx];
of << " " << arena.rhos_[fieldIdx];
}
Expand Down
15 changes: 8 additions & 7 deletions src/read_in_parameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,26 +355,27 @@ InitData read_in_parameters(std::string input_file) {
if (tempinput != "empty")
istringstream(tempinput) >> tempturn_on_QS;
parameter_list.turn_on_QS = tempturn_on_QS;
if (parameter_list.turn_on_QS == 1)
if (parameter_list.turn_on_QS == 1) {
parameter_list.alpha_max = 7;
else if (parameter_list.turn_on_rhob == 1)
} else if (parameter_list.turn_on_rhob == 1) {
parameter_list.alpha_max = 5;
else
} else {
parameter_list.alpha_max = 4;
}

int tempuse_BQ_ratios;
int tempuse_BQ_ratios = 0;
tempinput = Util::StringFind4(input_file, "use_BQ_ratios");
if (tempinput != "empty")
istringstream(tempinput) >> tempuse_BQ_ratios;
parameter_list.use_BQ_ratios = tempuse_BQ_ratios;

int tempuse_rhoQS_to_rhoB_ratios;
int tempuse_rhoQS_to_rhoB_ratios = 0;
tempinput = Util::StringFind4(input_file, "use_rhoQS_to_rhoB_ratios");
if (tempinput != "empty")
istringstream(tempinput) >> tempuse_rhoQS_to_rhoB_ratios;
parameter_list.use_rhoQS_to_rhoB_ratios = tempuse_rhoQS_to_rhoB_ratios;
double temp_ratio_q = 1;
double temp_ratio_s = 1;
double temp_ratio_q = 0.4;
double temp_ratio_s = 0.;
if(parameter_list.use_rhoQS_to_rhoB_ratios == 1){
tempinput = Util::StringFind4(input_file, "ratio_q");
if (tempinput != "empty")
Expand Down
10 changes: 5 additions & 5 deletions src/reconst.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@
Reconst::Reconst(const EOS &eosIn, const int echo_level_in, int beastMode) :
eos(eosIn),
max_iter(100),
LARGE(1e20),
v_critical(0.563624),
echo_level(echo_level_in) {
if (beastMode != 0) {
Expand Down Expand Up @@ -105,7 +104,7 @@ int Reconst::ReconstIt_velocity_Newton(ReconstCell &grid_p, double tau,
rhoq = J0Q/u[0];
rhos = J0S/u[0];

if (v_solution > v_critical) {
if (v_solution > v_critical && epsilon > 1e-6) {
// for large velocity, solve u0
double u0_solution = u[0];
int u0_status = solve_u0_Hybrid(u[0], T00, K00, M, J0B, J0Q, J0S,
Expand Down Expand Up @@ -474,7 +473,8 @@ void Reconst::reconst_velocity_f(const double v, const double T00,

double pressure = eos.get_pressure(epsilon, rhob, rhoq, rhos);

fv = v*(T00 + pressure) - M;
//fv = v*(T00 + pressure) - M;
fv = v - M/(T00 + pressure);
}


Expand Down Expand Up @@ -545,8 +545,8 @@ void Reconst::reconst_u0_f(const double u0, const double T00,

double pressure = eos.get_pressure(epsilon, rhob, rhoq, rhos);
const double temp1 = 1. - K00/((T00 + pressure)*(T00 + pressure));
//fu0 = u0 - 1./sqrt(temp1);
fu0 = u0*sqrt(temp1) - 1.;
fu0 = u0 - 1./sqrt(temp1);
//fu0 = u0*sqrt(temp1) - 1.;
}


Expand Down
4 changes: 1 addition & 3 deletions src/reconst.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,11 @@ class Reconst {
double rel_err;
double abs_err;

const double LARGE;

const double v_critical;
const int echo_level;

public:
Reconst() = default;
Reconst() = delete;
Reconst(const EOS &eos, const int echo_level_in, int beastMode);

ReconstCell ReconstIt_shell(double tau, const TJbVec &tauq_vec,
Expand Down
4 changes: 2 additions & 2 deletions src/transport_coeffs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@

using Util::hbarc;

TransportCoeffs::TransportCoeffs(const EOS &eosIn, const InitData &Data_in)
: DATA(Data_in), eos(eosIn),
TransportCoeffs::TransportCoeffs(const InitData &Data_in)
: DATA(Data_in),
shear_T_(DATA.T_dependent_shear_to_s),
shear_muB_(DATA.muB_dependent_shear_to_s),
bulk_T_(DATA.T_dependent_bulk_to_s) {
Expand Down
Loading

0 comments on commit c7587a1

Please sign in to comment.