SMASH Initial conditions - WIP

CMakeLists.txt

@@ -75,7 +75,7 @@ endif()
 
 string(APPEND CMAKE_CXX_FLAGS " -Wall")
 set(CMAKE_CXX_FLAGS_DEBUG "-g -O0")
-set(CMAKE_CXX_FLAGS_RELEASE "-g -O3")
+set(CMAKE_CXX_FLAGS_RELEASE "-g -O0")
 
 find_package(Git)
 if(Git_FOUND)

compile.sh

@@ -5,7 +5,7 @@ cd build
 rm -fr *
 cmake .. -DCMAKE_BUILD_TYPE=Release
 #CXX=g++-12 cmake .. -DCMAKE_BUILD_TYPE=Release
-#CXX=g++-12 cmake .. -DCMAKE_BUILD_TYPE=Debug
+#cmake .. -DCMAKE_BUILD_TYPE=Debug
 make -j4
 make install
 cd ..

input_params/RuRuSMASH200

@@ -0,0 +1,80 @@
+eho_level  1
+beastMode  0
+mode  2
+Initial_profile 31 
+ecm  200.0
+string_source_sigma_x  0.4
+string_source_sigma_eta  0.5
+stringTransverseShiftFrac  1.0
+stringPreEqFlowFactor  0.5
+Initial_Distribution_SMASH_filename initial/SMASH_IC.oscar
+event_id_SMASH_output 1
+average_SMASH_events 0
+extended_SMASH_output 1
+reject_SMASH_spectators 1
+s_factor  1.0
+Initial_time_tau_0  1.0
+Delta_Tau  0.02
+Total_evolution_time_tau  30.0
+boost_invariant  0
+Eta_grid_size  16.0
+Grid_size_in_eta  64
+X_grid_size_in_fm  30 
+Y_grid_size_in_fm  30 
+Grid_size_in_x  96 
+Grid_size_in_y  96
+gridPadding  3
+EOS_to_use  20 
+Shear_to_S_ratio 0.08 
+quest_revert_strength  1.0
+Viscosity_Flag_Yes_1_No_0  1
+Include_Shear_Visc_Yes_1_No_0 1
+shear_viscosity_3_eta_over_s_T_kink_in_GeV 0.17
+shear_viscosity_3_eta_over_s_low_T_slope_in_GeV -0.0
+shear_viscosity_3_eta_over_s_high_T_slope_in_GeV 0.0
+shear_viscosity_3_eta_over_s_at_kink 0.08
+muB_dependent_Shear_to_S_ratio 10
+shear_muBDep_alpha 0.7
+shear_muBDep_slope 2.0
+shear_muBDep_scale 0.6
+Include_Bulk_Visc_Yes_1_No_0  1
+T_dependent_zeta_over_s  10 
+bulk_viscosity_10_max 0.10  
+bulk_viscosity_10_T_peak 0.17
+bulk_viscosity_10_T_peak_muBcurv -0.15
+bulk_viscosity_10_width_high 0.08,
+bulk_viscosity_10_width_low 0.010,
+Include_second_order_terms  1
+Include_vorticity_terms  0
+Include_Rhob_Yes_1_No_0  1
+Include_QS_Yes_1_No_0  1
+use_BQ_ratios  0
+turn_on_baryon_diffusion  0
+kappa_coefficient  0.4
+output_hydro_debug_info  0
+output_evolution_data  2
+output_movie_flag  0
+output_evolution_T_cut  0.145
+outputBinaryEvolution  1
+output_evolution_every_N_eta  1
+output_evolution_every_N_x  1
+output_evolution_every_N_y  1
+output_evolution_every_N_timesteps  10
+Do_FreezeOut_Yes_1_No_0  1
+Do_FreezeOut_lowtemp  1
+freeze_out_tau_start_max  2
+freeze_out_method  4
+freeze_surface_in_binary  1
+average_surface_over_this_many_time_steps  10
+freeze_Ncell_x_step  1
+freeze_Ncell_eta_step  1
+freeze_eps_flag  0
+N_freeze_out  1
+eps_switch  0.35
+eps_freeze_max  0.5
+eps_freeze_min  0.5
+use_eps_for_freeze_out  1
+nucleon_width  0.4
+yL_frac  0.15
+output_vorticity  0
+EndOfData

src/CMakeLists.txt

@@ -9,6 +9,7 @@ set (SOURCES
     hydro_source_base.cpp
     hydro_source_strings.cpp
     hydro_source_ampt.cpp
+    hydro_source_smash.cpp
     hydro_source_TATB.cpp
     pretty_ostream.cpp
     grid_info.cpp

src/advance.cpp

@@ -154,7 +154,10 @@ void Advance::FirstRKStepT(
 
         if (DATA.turn_on_QS == 1) {
             qi_source[5] = (
-                tau_rk*hydro_source_terms_ptr->get_hydro_rhoQ_source(
+                tau_rk*hydro_source_terms_ptr->get_hydro_rhoq_source(
+                            tau_rk, x_local, y_local, eta_s_local, u_local));
+            qi_source[6] = (
+                tau_rk*hydro_source_terms_ptr->get_hydro_rhos_source(
                             tau_rk, x_local, y_local, eta_s_local, u_local));
         }
     }

src/data.h

@@ -35,6 +35,14 @@ typedef struct init_data {
     // parameters for Initial_profile == 11 || 111
     std::string initName_TA;  //!< filename for nuclear thickness function TA
     std::string initName_TB;  //!< filename for nuclear thickness function TB
+
+    // SMASH Oscar reading parameters
+    std::string initName_SMASH; //! filename for smash input file
+    int event_id_SMASH_output; //! selected event index. 
+    int average_SMASH_events; //! Average over events until reaching event_id_SMASH
+    int extended_SMASH_output; //! Consider extended Oscar input. 
+    int reject_SMASH_spectators; //! Do not consider spectators. Default: 1.
+
     std::string initName_participants;  //!< filename for participant nucleons
     double nucleonWidth;
 

src/hydro_source_base.cpp

@@ -40,7 +40,48 @@ double HydroSourceBase::get_hydro_rhob_source_before_tau(
     int n_tau_steps = static_cast<int>((tau - tau0)/dtau);
     for (int i = 0; i < n_tau_steps; i++) {
         const double tau_local = tau0 + (i + 0.5)*dtau;
-        const double res_local = get_hydro_rhob_source(
+        const double res_local = get_hydro_rhob_source(tau_local, x, y, eta_s, u);
+        res += tau_local*res_local*dtau;
+    }
+
+    return(res/tau);
+}
+
+double HydroSourceBase::get_hydro_rhoq_source_before_tau(
+                                const double tau, const double x,
+                                const double y, const double eta_s) const {
+    FlowVec u = {0};
+    u[0] = 1.0;
+
+    double res  = 0.;
+    double tau0 = 0.0;
+    double dtau = 0.005;
+
+    int n_tau_steps = static_cast<int>((tau - tau0)/dtau);
+    for (int i = 0; i < n_tau_steps; i++) {
+        const double tau_local = tau0 + (i + 0.5)*dtau;
+        const double res_local = get_hydro_rhoq_source(
+                                                tau_local, x, y, eta_s, u);
+        res += tau_local*res_local*dtau;
+    }
+
+    return(res/tau);
+}
+
+double HydroSourceBase::get_hydro_rhos_source_before_tau(
+                                const double tau, const double x,
+                                const double y, const double eta_s) const {
+    FlowVec u = {0};
+    u[0] = 1.0;
+
+    double res  = 0.;
+    double tau0 = 0.0;
+    double dtau = 0.005;
+
+    int n_tau_steps = static_cast<int>((tau - tau0)/dtau);
+    for (int i = 0; i < n_tau_steps; i++) {
+        const double tau_local = tau0 + (i + 0.5)*dtau;
+        const double res_local = get_hydro_rhos_source(
                                                 tau_local, x, y, eta_s, u);
         res += tau_local*res_local*dtau;
     }

src/hydro_source_base.h

@@ -58,7 +58,12 @@ class HydroSourceBase {
 
     //! this function returns the net electric charge density source term rho
     //! at a given point (tau, x, y, eta_s)
-    virtual double get_hydro_rhoQ_source(const double tau, const double x,
+    virtual double get_hydro_rhoq_source(const double tau, const double x,
+                                 const double y, const double eta_s,
+                                 const FlowVec &u_mu) const {
+        return(0.0);
+    }
+    virtual double get_hydro_rhos_source(const double tau, const double x,
                                  const double y, const double eta_s,
                                  const FlowVec &u_mu) const {
         return(0.0);
@@ -75,6 +80,12 @@ class HydroSourceBase {
     double get_hydro_rhob_source_before_tau(const double tau, const double x,
                                             const double y,
                                             const double eta_s) const;
+    double get_hydro_rhoq_source_before_tau(const double tau, const double x,
+                                            const double y,
+                                            const double eta_s) const;
+    double get_hydro_rhos_source_before_tau(const double tau, const double x,
+                                            const double y,
+                                            const double eta_s) const;
 
     virtual void prepare_list_for_current_tau_frame(const double tau_local) {}
 };