Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
causal_liquifier.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file causal_liquifier.cc
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2018
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 
16 #include "CausalLiquefier.h"
17 #include "LiquefierBase.h"
18 #include "gtest/gtest.h"
19 #include <iostream>
20 
21 using namespace Jetscape;
22 
23 // check coordinate transformation
24 TEST(CausalLiquifierTest, TEST_COORDINATES){
25 
26  CausalLiquefier lqf(0.3,0.3,0.3,0.3);
27 
28  // for the transformation from tau-eta to t-z (configuration space)
29  EXPECT_DOUBLE_EQ(0.0, lqf.get_t(0.0,0.5));
30  EXPECT_DOUBLE_EQ(5.0, lqf.get_t(5.0,0.0));
31  EXPECT_DOUBLE_EQ(0.0, lqf.get_z(0.0,0.5));
32  EXPECT_DOUBLE_EQ(0.0, lqf.get_z(5.0,0.0));
33 
34  for(double tau=0.1;tau<0.3;tau+=0.1){
35  for(double eta=0.0;eta<0.2;eta+=0.1){
36  double t=lqf.get_t(tau,eta);
37  double z=lqf.get_z(tau,eta);
38  EXPECT_DOUBLE_EQ(tau*tau, t*t-z*z);
39  EXPECT_DOUBLE_EQ(exp(2.0*eta), (t+z)/(t-z));
40  }
41  }
42 
43  // for the transformation from t-z to tau-eta (momentum)
44  EXPECT_DOUBLE_EQ(0.0, lqf.get_ptau(0.0,0.0,0.5));
45  EXPECT_DOUBLE_EQ(5.0, lqf.get_ptau(5.0,0.0,0.0));
46  EXPECT_DOUBLE_EQ(0.0, lqf.get_peta(0.0,0.0,0.5));
47  EXPECT_DOUBLE_EQ(5.0, lqf.get_peta(0.0,5.0,0.0));
48 
49 }
50 
51 // check causality
52 TEST(CausalLiquifierTest, TEST_CAUSALITY){
53 
54  CausalLiquefier lqf(0.3,0.3,0.3,0.3);
55 
56  // check zero source outside of the cousal area
57  double time = lqf.tau_delay;
58  double r_bound = lqf.c_diff * time;
59  for(int scale = 1.0; scale < 5.0; scale++ ){
60  double r_out = r_bound * scale;
61  EXPECT_DOUBLE_EQ(0.0, lqf.rho_smooth(time,r_out));
62  EXPECT_DOUBLE_EQ(0.0, lqf.rho_delta(time,r_out));
63  EXPECT_DOUBLE_EQ(0.0, lqf.kernel_rho(time,r_out));
64  }
65 
66  // check zero source before the deposition time
67  time = 0.99*(lqf.tau_delay-0.5*lqf.dtau);
68  std::array<Jetscape::real, 4> jmu = {0.0,0.0,0.0,0.0};
69  std::array<Jetscape::real, 4> droplet_xmu = {0.0, 0.0, 0.0, 0.0};
70  std::array<Jetscape::real, 4> droplet_pmu = {1.0, 1.0, 0.0, 0.0};
71  Droplet drop_i(droplet_xmu, droplet_pmu);
72  lqf.smearing_kernel( time, 0.0, 0.0, 0.0, drop_i, jmu);
73  EXPECT_DOUBLE_EQ(0.0, jmu[0]);
74 
75  // check zero source after the deposition time
76  time = 1.01*(lqf.tau_delay+0.5*lqf.dtau);
77  lqf.smearing_kernel( time, 0.0, 0.0, 0.0, drop_i, jmu);
78  EXPECT_DOUBLE_EQ(0.0, jmu[0]);
79 }
80 
81 
82 
83 // check conservation law
84 TEST(CausalLiquifierTest, TEST_CONSERVATION){
85 
86  double dr = 0.005; //in [fm]
87  double dt = 0.005; //in [fm]
88 
89  CausalLiquefier lqf(dt,dr,dr,dr);
90 
91  for(double t=0.3; t < 3.0; t+=dt){
92  double integrated_value = 0.0;
93  for(double r = 0.5*dr; r<1.5*t; r+=dr){
94  integrated_value
95  += 4.0*M_PI*r*r*dr*lqf.kernel_rho(t,r);
96  }
97  // require conservation within 5%
98  EXPECT_NEAR(1.0,integrated_value,0.05);
99  }
100 
101 }
102 
103 // check conservation law
104 TEST(CausalLiquifierTest, TEST_GRID_CARTESIAN_CONSERVATION){
105 
106  double ll[3] = {0.05,0.1,0.3};
107 
108  for( int i=0; i<0; i++){
109 
110  double dt = 0.1;
111  double l = ll[i];
112  double dx = l;
113  double dy = l;
114  double dz = l;
115  int n_cells = 5.0/l;
116  double dV = dx*dx*dz;
117 
118  CausalLiquefier lqf(dt,dx,dy,dz);
119 
120  //std::o/Users/yasukitachibana/Dropbox (Personal)/Codes/Release2.2CandidateLiquefierUpdate/examples/unittests/causal_liquifier.ccfstream ofs;
121  //string filename = "liq_cons_test_dx" + std::to_string(int(l*1000)) + "_tr3000_d2400_delta300.txt";
122  //ofs.open(filename.c_str(), std::ios_base::out);
123 
124 
125  for(double t=0.5; t < 3.0; t+=dt){
126  double integrated_value = 0.0;
127  for(int ix=0; ix<n_cells; ix++ ){
128  double x = ( ix - 0.5*double(n_cells) )*dx;
129  for(int iy=0; iy<n_cells; iy++ ){
130  double y = ( iy - 0.5*double(n_cells) )*dy;
131  for(int iz=0; iz<n_cells; iz++ ){
132  double z = ( iz - 0.5*double(n_cells) )*dz;
133 
134  double r = sqrt(x*x+y*y+z*z);
135  integrated_value += lqf.kernel_rho(t,r)*dV;
136 
137 
138  }
139  }
140  }
141  // require conservation within 5%
142  //ofs << t << " " << integrated_value <<"\n";
143  EXPECT_NEAR(1.0,integrated_value,0.05);
144  //JSINFO <<"t="<<t <<" " <<integrated_value;
145  }
146  //ofs.close();
147  }
148 }
149 
150 
151 // check conservation law
152 TEST(CausalLiquifierTest, TEST_GRID_TAU_ETA_CONSERVATION){
153 
154 
155  std::array<Jetscape::real, 4> x_in = {1.0, 0.0, 0.0, 1.0};// in tau-eta
156  std::array<Jetscape::real, 4> p_in = {1.0, 1.0, 1.0, 1.0};// in Cartesian
157  Droplet a_drop(x_in,p_in);
158 
159 
160  double dtau = 0.1;
161  double dl = 0.05;
162  double dx = dl;
163  double dy = dl;
164  double deta = 0.05;
165  int n_xy = 8.0/dl;
166  int n_eta = 5.0/deta;
167 
168  CausalLiquefier lqf(dtau,dx,dy,deta);
169 
170 
171 // std::ofstream ofs;
172 // string filename = "liq_cons_tau_eta_test_tdep"+std::to_string(int(x_in[0]))+"_etadep"+std::to_string(int(x_in[3]))+"_dxy" + std::to_string(int(dl*1000)) +
173 // + "_deta" + std::to_string(int(deta*1000)) + "_tr100_d80_delta100.txt";
174 //
175 // ofs.open(filename.c_str(), std::ios_base::out);
176 
177  double dt = 0.5;
178  for(double t=0.2; t < 4.0; t+=dt){
179  double integrated_value = 0.0;
180  double tau_delay = t;
181  lqf.set_t_delay(tau_delay);
182  // in tau-eta
183  std::array<Jetscape::real, 4> total_pmu = {0.0, 0.0, 0.0, 0.0};
184  std::array<Jetscape::real, 4> x_hydro = {0.0, 0.0, 0.0, 0.0};
185  x_hydro[0] = x_in[0]+tau_delay;
186  double dvolume = x_hydro[0]*dx*dy*deta;
187 
188  for(int ix=0; ix<n_xy; ix++ ){
189  x_hydro[1] = ( ix - 0.5*double(n_xy) )*dx;
190  for(int iy=0; iy<n_xy; iy++ ){
191  x_hydro[2] = ( iy - 0.5*double(n_xy) )*dy;
192  for(int ieta=0; ieta<n_eta; ieta++ ){
193  x_hydro[3] = ( ieta - 0.5*double(n_eta) )*deta;
194 
195 
196  std::array<Jetscape::real, 4> jmu = {0.0, 0.0, 0.0, 0.0};
197 
198  lqf.smearing_kernel(x_hydro[0],x_hydro[1],x_hydro[2],x_hydro[3],a_drop, jmu);
199 
200  total_pmu[0] += dtau*(jmu[0]*cosh(x_hydro[3]) + jmu[3]*sinh(x_hydro[3]))*dvolume;
201  total_pmu[1] += dtau*jmu[1]*dvolume;
202  total_pmu[2] += dtau*jmu[2]*dvolume;
203  total_pmu[3] += dtau*(jmu[0]*sinh(x_hydro[3]) + jmu[3]*cosh(x_hydro[3]))*dvolume;
204 
205 
206  integrated_value += dtau*jmu[1]*dvolume;
207 
208  }
209  }
210  }
211 // ofs << t << " " << integrated_value <<"\n";
212 // JSINFO
213 // << total_pmu[0] << " "
214 // << total_pmu[1] << " "
215 // << total_pmu[2] << " "
216 // << total_pmu[3];
217 // ofs.close();
218  if(t>3.){
219  EXPECT_NEAR(1.0,integrated_value,0.05);
220  }
221  }
222 }