examples/unittests/causal_liquifier.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 #include "CausalLiquefier.h"
0017 #include "LiquefierBase.h"
0018 #include "gtest/gtest.h"
0019 #include <iostream>
0020
0021 using namespace Jetscape;
0022
0023 // check coordinate transformation
0024 TEST(CausalLiquifierTest, TEST_COORDINATES){
0025
0026     CausalLiquefier lqf(0.3,0.3,0.3,0.3);
0027
0028     // for the transformation from tau-eta to t-z (configuration space)
0029     EXPECT_DOUBLE_EQ(0.0, lqf.get_t(0.0,0.5));
0030     EXPECT_DOUBLE_EQ(5.0, lqf.get_t(5.0,0.0));
0031     EXPECT_DOUBLE_EQ(0.0, lqf.get_z(0.0,0.5));
0032     EXPECT_DOUBLE_EQ(0.0, lqf.get_z(5.0,0.0));
0033
0034     for(double tau=0.1;tau<0.3;tau+=0.1){
0035         for(double eta=0.0;eta<0.2;eta+=0.1){
0036             double t=lqf.get_t(tau,eta);
0037             double z=lqf.get_z(tau,eta);
0038             EXPECT_DOUBLE_EQ(tau*tau, t*t-z*z);
0039             EXPECT_DOUBLE_EQ(exp(2.0*eta), (t+z)/(t-z));
0040         }
0041     }
0042
0043     // for the transformation from t-z to tau-eta (momentum)
0044     EXPECT_DOUBLE_EQ(0.0, lqf.get_ptau(0.0,0.0,0.5));
0045     EXPECT_DOUBLE_EQ(5.0, lqf.get_ptau(5.0,0.0,0.0));
0046     EXPECT_DOUBLE_EQ(0.0, lqf.get_peta(0.0,0.0,0.5));
0047     EXPECT_DOUBLE_EQ(5.0, lqf.get_peta(0.0,5.0,0.0));
0048
0049 }
0050
0051 // check causality
0052 TEST(CausalLiquifierTest, TEST_CAUSALITY){
0053
0054     CausalLiquefier lqf(0.3,0.3,0.3,0.3);
0055
0056     // check zero source outside of the cousal area
0057     double time = lqf.tau_delay;
0058     double r_bound = lqf.c_diff * time;
0059     for(int scale = 1.0; scale < 5.0; scale++ ){
0060         double r_out = r_bound * scale;
0061         EXPECT_DOUBLE_EQ(0.0, lqf.rho_smooth(time,r_out));
0062         EXPECT_DOUBLE_EQ(0.0, lqf.rho_delta(time,r_out));
0063         EXPECT_DOUBLE_EQ(0.0, lqf.kernel_rho(time,r_out));
0064     }
0065
0066     // check zero source before the deposition time
0067     time = 0.99*(lqf.tau_delay-0.5*lqf.dtau);
0068     std::array<Jetscape::real, 4> jmu = {0.0,0.0,0.0,0.0};
0069     std::array<Jetscape::real, 4> droplet_xmu = {0.0, 0.0, 0.0, 0.0};
0070     std::array<Jetscape::real, 4> droplet_pmu = {1.0, 1.0, 0.0, 0.0};
0071     Droplet drop_i(droplet_xmu, droplet_pmu);
0072     lqf.smearing_kernel( time, 0.0, 0.0, 0.0, drop_i, jmu);
0073     EXPECT_DOUBLE_EQ(0.0, jmu[0]);
0074
0075     // check zero source after the deposition time
0076     time = 1.01*(lqf.tau_delay+0.5*lqf.dtau);
0077     lqf.smearing_kernel( time, 0.0, 0.0, 0.0, drop_i, jmu);
0078     EXPECT_DOUBLE_EQ(0.0, jmu[0]);
0079 }
0080
0081
0082
0083 // check conservation law
0084 TEST(CausalLiquifierTest, TEST_CONSERVATION){
0085
0086     double dr = 0.005; //in [fm]
0087     double dt = 0.005; //in [fm]
0088
0089     CausalLiquefier lqf(dt,dr,dr,dr);
0090
0091     for(double t=0.3; t < 3.0; t+=dt){
0092         double integrated_value = 0.0;
0093         for(double r = 0.5*dr; r<1.5*t; r+=dr){
0094             integrated_value
0095             += 4.0*M_PI*r*r*dr*lqf.kernel_rho(t,r);
0096         }
0097         // require conservation within 5%
0098         EXPECT_NEAR(1.0,integrated_value,0.05);
0099     }
0100
0101 }
0102
0103 // check conservation law
0104 TEST(CausalLiquifierTest, TEST_GRID_CARTESIAN_CONSERVATION){
0105
0106     double ll[3] = {0.05,0.1,0.3};
0107
0108     for( int i=0; i<0; i++){
0109
0110         double dt = 0.1;
0111         double l = ll[i];
0112         double dx = l;
0113         double dy = l;
0114         double dz = l;
0115         int n_cells = 5.0/l;
0116         double dV = dx*dx*dz;
0117
0118         CausalLiquefier lqf(dt,dx,dy,dz);
0119
0120         //std::o/Users/yasukitachibana/Dropbox (Personal)/Codes/Release2.2CandidateLiquefierUpdate/examples/unittests/causal_liquifier.ccfstream ofs;
0121         //string filename = "liq_cons_test_dx" + std::to_string(int(l*1000)) + "_tr3000_d2400_delta300.txt";
0122         //ofs.open(filename.c_str(), std::ios_base::out);
0123
0124
0125         for(double t=0.5; t < 3.0; t+=dt){
0126             double integrated_value = 0.0;
0127             for(int ix=0; ix<n_cells; ix++ ){
0128                 double x = ( ix - 0.5*double(n_cells) )*dx;
0129                 for(int iy=0; iy<n_cells; iy++ ){
0130                     double y = ( iy - 0.5*double(n_cells) )*dy;
0131                     for(int iz=0; iz<n_cells; iz++ ){
0132                         double z = ( iz - 0.5*double(n_cells) )*dz;
0133
0134                         double r = sqrt(x*x+y*y+z*z);
0135                         integrated_value += lqf.kernel_rho(t,r)*dV;
0136
0137
0138                     }
0139                 }
0140             }
0141             // require conservation within 5%
0142             //ofs << t << " " << integrated_value <<"\n";
0143             EXPECT_NEAR(1.0,integrated_value,0.05);
0144             //JSINFO <<"t="<<t <<" " <<integrated_value;
0145         }
0146         //ofs.close();
0147     }
0148 }
0149
0150
0151 // check conservation law
0152 TEST(CausalLiquifierTest, TEST_GRID_TAU_ETA_CONSERVATION){
0153
0154
0155     std::array<Jetscape::real, 4> x_in = {1.0, 0.0, 0.0, 1.0};// in tau-eta
0156     std::array<Jetscape::real, 4> p_in = {1.0, 1.0, 1.0, 1.0};// in Cartesian
0157     Droplet a_drop(x_in,p_in);
0158
0159
0160     double dtau = 0.1;
0161     double dl = 0.05;
0162     double dx = dl;
0163     double dy = dl;
0164     double deta = 0.05;
0165     int n_xy = 8.0/dl;
0166     int n_eta = 5.0/deta;
0167
0168     CausalLiquefier lqf(dtau,dx,dy,deta);
0169
0170
0171 //    std::ofstream ofs;
0172 //    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)) +
0173 //    + "_deta" + std::to_string(int(deta*1000)) + "_tr100_d80_delta100.txt";
0174 //
0175 //    ofs.open(filename.c_str(), std::ios_base::out);
0176
0177     double dt = 0.5;
0178     for(double t=0.2; t < 4.0; t+=dt){
0179         double integrated_value = 0.0;
0180         double tau_delay = t;
0181         lqf.set_t_delay(tau_delay);
0182         // in tau-eta
0183         std::array<Jetscape::real, 4> total_pmu = {0.0, 0.0, 0.0, 0.0};
0184         std::array<Jetscape::real, 4> x_hydro = {0.0, 0.0, 0.0, 0.0};
0185         x_hydro[0] = x_in[0]+tau_delay;
0186         double dvolume = x_hydro[0]*dx*dy*deta;
0187
0188         for(int ix=0; ix<n_xy; ix++ ){
0189             x_hydro[1] = ( ix - 0.5*double(n_xy) )*dx;
0190             for(int iy=0; iy<n_xy; iy++ ){
0191                 x_hydro[2] = ( iy - 0.5*double(n_xy) )*dy;
0192                 for(int ieta=0; ieta<n_eta; ieta++ ){
0193                     x_hydro[3] = ( ieta - 0.5*double(n_eta) )*deta;
0194
0195
0196                     std::array<Jetscape::real, 4> jmu = {0.0, 0.0, 0.0, 0.0};
0197
0198                     lqf.smearing_kernel(x_hydro[0],x_hydro[1],x_hydro[2],x_hydro[3],a_drop, jmu);
0199
0200                     total_pmu[0] += dtau*(jmu[0]*cosh(x_hydro[3]) + jmu[3]*sinh(x_hydro[3]))*dvolume;
0201                     total_pmu[1] += dtau*jmu[1]*dvolume;
0202                     total_pmu[2] += dtau*jmu[2]*dvolume;
0203                     total_pmu[3] += dtau*(jmu[0]*sinh(x_hydro[3]) + jmu[3]*cosh(x_hydro[3]))*dvolume;
0204
0205
0206                     integrated_value += dtau*jmu[1]*dvolume;
0207
0208                 }
0209             }
0210         }
0211 //        ofs << t << " " << integrated_value <<"\n";
0212 //        JSINFO
0213 //        << total_pmu[0] << " "
0214 //        << total_pmu[1] << " "
0215 //        << total_pmu[2] << " "
0216 //        << total_pmu[3];
0217 //        ofs.close();
0218         if(t>3.){
0219             EXPECT_NEAR(1.0,integrated_value,0.05);
0220         }
0221     }
0222 }