src/liquefier/CausalLiquefier.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015 // -----------------------------------------
0016 // This is a causal liquefier with the JETSCAPE framework
0017 // -----------------------------------------
0018
0019 #include "CausalLiquefier.h"
0020 #include "JetScapeLogger.h"
0021 #include "JetScapeXML.h"
0022 #include <cfloat>
0023
0024 namespace Jetscape {
0025
0026
0027 CausalLiquefier::CausalLiquefier(){
0028     VERBOSE(8);
0029     dtau = 0.6;
0030     dx = 0.3;
0031     dy = 0.3;
0032     deta = 0.3;
0033     tau_delay = 2.0;
0034     time_relax = 0.1;
0035     d_diff = 0.08;
0036     width_delta = 0.1;
0037     Init();// Get values of parameters from XML
0038     c_diff = sqrt(d_diff/time_relax);
0039     gamma_relax = 0.5/time_relax;
0040     if( c_diff > 1.0 ){
0041         JSWARN << "Bad Signal Velocity in CausalLiquefier";
0042     }
0043 //    else{
0044 //        //for debug
0045 //        JSINFO << "c_diff = " << c_diff;
0046 //    }
0047
0048 }
0049
0050 CausalLiquefier::CausalLiquefier(double dtau_in, double dx_in, double dy_in, double deta_in){
0051     VERBOSE(8);
0052     JSINFO<<"Initialize CausalLiquefier (for Unit Test) ...";
0053     dtau = dtau_in;
0054     dx = dx_in;
0055     dy = dy_in;
0056     deta = deta_in;
0057     tau_delay = 2.0;
0058     time_relax = 0.1;
0059     d_diff = 0.08;
0060     width_delta = 0.1;
0061
0062     c_diff = sqrt(d_diff/time_relax);
0063     gamma_relax = 0.5/time_relax;
0064
0065     JSINFO
0066     << "<CausalLiquefier> Fluid Time Step and Cell Size: dtau="
0067     << dtau << " fm, dx="
0068     << dx << " fm, dy="
0069     << dy << " fm, deta="
0070     << deta;
0071     JSINFO
0072     << "<CausalLiquefier> Parameters: tau_delay="
0073     << tau_delay << " fm, time_relax="
0074     << time_relax << " fm, d_diff="
0075     << d_diff << " fm, width_delta="
0076     << width_delta <<" fm";
0077 }
0078
0079 void CausalLiquefier::Init(){
0080
0081     // Initialize parameter with values in XML
0082     JSINFO<<"Initialize CausalLiquefier ...";
0083
0084     dtau = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "dtau"});
0085     dx = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "dx"});
0086     dy = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "dy"});
0087     deta = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "deta"});
0088     tau_delay = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "tau_delay"});// in [fm]
0089     time_relax = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "time_relax"});// in [fm]
0090     d_diff = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "d_diff"});// in [fm]
0091     width_delta = JetScapeXML::Instance()->GetElementDouble({"Liquefier", "CausalLiquefier", "width_delta"});// in [fm]
0092
0093     // for debug
0094 //    JSINFO
0095 //    << "<CausalLiquefier> Fluid Time Step and Cell Size: dtau="
0096 //    << dtau << " fm, dx="
0097 //    << dx << " fm, dy="
0098 //    << dy << " fm, deta="
0099 //    << deta;
0100 //    JSINFO
0101 //    << "<CausalLiquefier> Parameters: tau_delay="
0102 //    << tau_delay << " fm, time_relax="
0103 //    << time_relax << " fm, d_diff="
0104 //    << d_diff << " fm, width_delta="
0105 //    << width_delta <<" fm";
0106
0107 }
0108
0109 //CausalLiquefier::~CausalLiquefier(){};
0110
0111
0112 void CausalLiquefier::smearing_kernel(
0113         Jetscape::real tau, Jetscape::real x, Jetscape::real y,
0114         Jetscape::real eta, const Droplet drop_i,
0115         std::array<Jetscape::real, 4> &jmu) const {
0116
0117     jmu = {0., 0, 0, 0};//source in tau-eta coordinates
0118
0119     const auto p_drop = drop_i.get_pmu();
0120     auto x_drop = drop_i.get_xmu();// position of the doloplet in the Cartesian coodinates
0121     double tau_drop = x_drop[0];
0122     double eta_drop = x_drop[3];
0123     x_drop[0] = get_t(tau_drop, eta_drop);
0124     x_drop[3] = get_z(tau_drop, eta_drop);
0125
0126     if( tau - 0.5*dtau <= tau_drop + tau_delay &&
0127        tau + 0.5*dtau > tau_drop + tau_delay ){
0128
0129         double t = get_t(tau, eta);// position in the fluid in the Cartesian coodinates
0130         double z = get_z(tau, eta);
0131
0132         double delta_t = t - x_drop[0];
0133         double delta_r = sqrt((x-x_drop[1])*(x-x_drop[1])+(y-x_drop[2])*(y-x_drop[2])+(z-x_drop[3])*(z-x_drop[3]));
0134
0135         // get solutions of the diffusion equation in the Cartesian coordinates
0136         double jt = kernel_rho( delta_t, delta_r )/dtau;
0137         double jz;
0138         if( delta_r <= DBL_MIN ){
0139             jz = 0.0;
0140         }else{
0141             jz = ((z-x_drop[3])/delta_r)*kernel_j( delta_t, delta_r )/dtau;
0142         }
0143         // get flux for the constant-tau surface
0144         double jtau = get_ptau(jt, jz, eta);
0145
0146         // get source in tau-eta coordinates
0147         jmu[0] = jtau*get_ptau(p_drop[0], p_drop[3], eta);
0148         jmu[1] = jtau*p_drop[1];
0149         jmu[2] = jtau*p_drop[2];
0150         jmu[3] = jtau*get_peta(p_drop[0], p_drop[3], eta);
0151
0152     }
0153
0154 }
0155
0156
0157 //Charge density rho in causal diffusion
0158 double CausalLiquefier::kernel_rho(double t, double r) const {
0159     return dumping(t)*(rho_smooth(t, r)+rho_delta(t, r));
0160 }
0161
0162 //Radial component of current j in causal diffusion
0163 double CausalLiquefier::kernel_j(double t, double r) const {
0164     return dumping(t)*(j_smooth(t, r)+j_delta(t, r));
0165 }
0166
0167 //Dumping factor in solutions of rho and j
0168 double CausalLiquefier::dumping(double t) const {
0169     return exp(-gamma_relax*t)/(4.0*M_PI);
0170 }
0171
0172 //Smooth component of rho
0173 double CausalLiquefier::rho_smooth(double t, double r) const {
0174     if( r < c_diff*t ){
0175         double u = sqrt( c_diff*c_diff*t*t - r*r );
0176         double x = gamma_relax*u/c_diff; // unitless
0177
0178         double i1 = gsl_sf_bessel_I1(x)/(c_diff*u);
0179         double i2 = gsl_sf_bessel_In(2,x)*t/u/u;
0180         double f = gamma_relax*gamma_relax/c_diff;
0181
0182         return f * (i1+i2);
0183     }else{
0184         return 0.0;
0185     }
0186 }
0187
0188 //Smooth component of j
0189 double CausalLiquefier::j_smooth(double t, double r) const {
0190     if( r < c_diff*t ){
0191         double u = sqrt( c_diff*c_diff*t*t - r*r );
0192         double x = gamma_relax*u/c_diff; // unitless
0193
0194         double i2 = gsl_sf_bessel_In(2,x)*r/u/u;
0195         double f = gamma_relax*gamma_relax/c_diff;
0196
0197         return f * i2;
0198
0199     }else{
0200         return 0.0;
0201     }
0202 }
0203
0204 //Wave front component of rho
0205 double CausalLiquefier::rho_delta(double t, double r) const {
0206
0207     double r_w = width_delta;
0208     if( c_diff*t <= width_delta ){
0209         r_w = c_diff*t;
0210     }
0211
0212     if( r >= c_diff*t - r_w && r < c_diff*t ){
0213         double x = gamma_relax*t;
0214         return (1.0 + x + x*x/2.0)/r_w/r/r;
0215     }else{
0216         return 0.0;
0217     }
0218
0219 }
0220
0221 //Wave front component of j
0222 double CausalLiquefier::j_delta(double t, double r) const {
0223     return c_diff*rho_delta(t, r);
0224 }
0225
0226 //Get Cartesian time t from tau and eta
0227 double CausalLiquefier::get_t(double tau, double eta)const{
0228     return tau*cosh(eta);
0229 }
0230
0231 //Get Cartesian coordinate z from tau and eta
0232 double CausalLiquefier::get_z(double tau, double eta)const{
0233         return tau*sinh(eta);
0234 }
0235
0236 //Lorentz Transformation to get tau component of four vector
0237 double CausalLiquefier::get_ptau(double p0, double p3, double eta)const{
0238         return p0*cosh(eta) - p3*sinh(eta);
0239 }
0240
0241 //Lorentz Transformation to get eta component of four vector
0242 double CausalLiquefier::get_peta(double p0, double p3, double eta)const{
0243         return p3*cosh(eta) - p0*sinh(eta);
0244 }
0245
0246 //For debug, Change tau_delay
0247 void CausalLiquefier::set_t_delay(double new_tau_delay){
0248         tau_delay = new_tau_delay;
0249 }
0250
0251 };