src/framework/SurfaceFinder.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 // This is a general basic class for a hyper-surface finder
0016
0017 #include <cmath>
0018 #include "RealType.h"
0019 #include "SurfaceFinder.h"
0020 #include "cornelius.h"
0021 #include "FluidEvolutionHistory.h"
0022 #include "JetScapeLogger.h"
0023
0024 namespace Jetscape {
0025
0026 SurfaceFinder::SurfaceFinder(const Jetscape::real T_in,
0027                              const EvolutionHistory &bulk_data)
0028     : bulk_info(bulk_data) {
0029
0030   T_cut = T_in;
0031   JSINFO << "Find a surface with temperature T = " << T_cut;
0032   boost_invariant = bulk_info.is_boost_invariant();
0033   if (boost_invariant) {
0034     JSINFO << "Hydro medium is boost invariant.";
0035   } else {
0036     JSINFO << "Hydro medium is not boost invariant.";
0037   }
0038   JSINFO << "Number of fluid cells = " << bulk_info.get_data_size();
0039 }
0040
0041 SurfaceFinder::~SurfaceFinder() { surface_cell_list.clear(); }
0042
0043 void SurfaceFinder::Find_full_hypersurface() {
0044   if (boost_invariant) {
0045     JSINFO << "Finding a 2+1D hyper-surface at T = " << T_cut << " GeV ...";
0046     Find_full_hypersurface_3D();
0047   } else {
0048     JSINFO << "Finding a 3+1D hyper-surface at T = " << T_cut << " GeV ...";
0049     Find_full_hypersurface_4D();
0050   }
0051 }
0052
0053 bool SurfaceFinder::check_intersect_3D(Jetscape::real tau, Jetscape::real x,
0054                                        Jetscape::real y, Jetscape::real dt,
0055                                        Jetscape::real dx, Jetscape::real dy,
0056                                        double ***cube) {
0057   bool intersect = true;
0058
0059   auto tau_low = tau - dt / 2.;
0060   auto tau_high = tau + dt / 2.;
0061   auto x_left = x - dx / 2.;
0062   auto x_right = x + dx / 2.;
0063   auto y_left = y - dy / 2.;
0064   auto y_right = y + dy / 2.;
0065
0066   auto fluid_cell = bulk_info.get(tau_low, x_left, y_left, 0.0);
0067   cube[0][0][0] = fluid_cell.temperature;
0068   fluid_cell = bulk_info.get(tau_low, x_left, y_right, 0.0);
0069   cube[0][0][1] = fluid_cell.temperature;
0070   fluid_cell = bulk_info.get(tau_low, x_right, y_left, 0.0);
0071   cube[0][1][0] = fluid_cell.temperature;
0072   fluid_cell = bulk_info.get(tau_low, x_right, y_right, 0.0);
0073   cube[0][1][1] = fluid_cell.temperature;
0074   fluid_cell = bulk_info.get(tau_high, x_left, y_left, 0.0);
0075   cube[1][0][0] = fluid_cell.temperature;
0076   fluid_cell = bulk_info.get(tau_high, x_left, y_right, 0.0);
0077   cube[1][0][1] = fluid_cell.temperature;
0078   fluid_cell = bulk_info.get(tau_high, x_right, y_left, 0.0);
0079   cube[1][1][0] = fluid_cell.temperature;
0080   fluid_cell = bulk_info.get(tau_high, x_right, y_right, 0.0);
0081   cube[1][1][1] = fluid_cell.temperature;
0082
0083   if ((T_cut - cube[0][0][0]) * (cube[1][1][1] - T_cut) < 0.0)
0084     if ((T_cut - cube[0][1][0]) * (cube[1][0][1] - T_cut) < 0.0)
0085       if ((T_cut - cube[0][1][1]) * (cube[1][0][0] - T_cut) < 0.0)
0086         if ((T_cut - cube[0][0][1]) * (cube[1][1][0] - T_cut) < 0.0)
0087           intersect = false;
0088
0089   return (intersect);
0090 }
0091
0092 void SurfaceFinder::Find_full_hypersurface_3D() {
0093   auto grid_tau0 = bulk_info.Tau0();
0094   auto grid_tauf = bulk_info.TauMax();
0095   auto grid_x0 = bulk_info.XMin();
0096   auto grid_y0 = bulk_info.YMin();
0097   ;
0098
0099   Jetscape::real grid_dt = 0.1;
0100   Jetscape::real grid_dx = 0.2;
0101   Jetscape::real grid_dy = 0.2;
0102
0103   const int dim = 3;
0104   double lattice_spacing[dim];
0105   lattice_spacing[0] = grid_dt;
0106   lattice_spacing[1] = grid_dx;
0107   lattice_spacing[2] = grid_dy;
0108
0109   std::unique_ptr<Cornelius> cornelius_ptr(new Cornelius());
0110   cornelius_ptr->init(dim, T_cut, lattice_spacing);
0111
0112   const int ntime = static_cast<int>((grid_tauf - grid_tau0) / grid_dt);
0113   const int nx = static_cast<int>(std::abs(2. * grid_x0) / grid_dx);
0114   const int ny = static_cast<int>(std::abs(2. * grid_y0) / grid_dy);
0115
0116   double ***cube = new double **[2];
0117   for (int i = 0; i < 2; i++) {
0118     cube[i] = new double *[2];
0119     for (int j = 0; j < 2; j++) {
0120       cube[i][j] = new double[2];
0121       for (int k = 0; k < 2; k++)
0122         cube[i][j][k] = 0.0;
0123     }
0124   }
0125
0126   for (int itime = 0; itime < ntime; itime++) {
0127     // loop over time evolution
0128     auto tau_local = grid_tau0 + (itime + 0.5) * grid_dt;
0129     for (int i = 0; i < nx; i++) {
0130       // loops over the transverse plane
0131       auto x_local = grid_x0 + (i + 0.5) * grid_dx;
0132       for (int j = 0; j < ny; j++) {
0133         auto y_local = grid_y0 + (j + 0.5) * grid_dy;
0134         bool intersect = check_intersect_3D(tau_local, x_local, y_local,
0135                                             grid_dt, grid_dx, grid_dy, cube);
0136         if (intersect) {
0137           cornelius_ptr->find_surface_3d(cube);
0138           for (int isurf = 0; isurf < cornelius_ptr->get_Nelements(); isurf++) {
0139             auto tau_center = (cornelius_ptr->get_centroid_elem(isurf, 0) +
0140                                tau_local - grid_dt / 2.);
0141             auto x_center = (cornelius_ptr->get_centroid_elem(isurf, 1) +
0142                              x_local - grid_dx / 2.);
0143             auto y_center = (cornelius_ptr->get_centroid_elem(isurf, 2) +
0144                              y_local - grid_dy / 2.);
0145
0146             auto da_tau = cornelius_ptr->get_normal_elem(isurf, 0);
0147             auto da_x = cornelius_ptr->get_normal_elem(isurf, 1);
0148             auto da_y = cornelius_ptr->get_normal_elem(isurf, 2);
0149
0150             auto fluid_cell =
0151                 bulk_info.get(tau_center, x_center, y_center, 0.0);
0152             auto surface_cell =
0153                 PrepareASurfaceCell(tau_center, x_center, y_center, 0.0, da_tau,
0154                                     da_x, da_y, 0.0, fluid_cell);
0155             surface_cell_list.push_back(surface_cell);
0156           }
0157         }
0158       }
0159     }
0160   }
0161
0162   for (int i = 0; i < 2; i++) {
0163     for (int j = 0; j < 2; j++)
0164       delete[] cube[i][j];
0165     delete[] cube[i];
0166   }
0167   delete[] cube;
0168 }
0169
0170 bool SurfaceFinder::check_intersect_4D(Jetscape::real tau, Jetscape::real x,
0171                                        Jetscape::real y, Jetscape::real eta,
0172                                        Jetscape::real dt, Jetscape::real dx,
0173                                        Jetscape::real dy, Jetscape::real deta,
0174                                        double ****cube) {
0175
0176   bool intersect = true;
0177
0178   auto tau_low = tau - dt / 2.;
0179   auto tau_high = tau + dt / 2.;
0180   auto x_left = x - dx / 2.;
0181   auto x_right = x + dx / 2.;
0182   auto y_left = y - dy / 2.;
0183   auto y_right = y + dy / 2.;
0184   auto eta_left = eta - deta / 2.;
0185   auto eta_right = eta + deta / 2.;
0186
0187   auto fluid_cell = bulk_info.get(tau_low, x_left, y_left, eta_left);
0188   cube[0][0][0][0] = fluid_cell.temperature;
0189   fluid_cell = bulk_info.get(tau_low, x_left, y_left, eta_right);
0190   cube[0][0][0][1] = fluid_cell.temperature;
0191   fluid_cell = bulk_info.get(tau_low, x_left, y_right, eta_left);
0192   cube[0][0][1][0] = fluid_cell.temperature;
0193   fluid_cell = bulk_info.get(tau_low, x_left, y_right, eta_right);
0194   cube[0][0][1][1] = fluid_cell.temperature;
0195   fluid_cell = bulk_info.get(tau_low, x_right, y_left, eta_left);
0196   cube[0][1][0][0] = fluid_cell.temperature;
0197   fluid_cell = bulk_info.get(tau_low, x_right, y_left, eta_right);
0198   cube[0][1][0][1] = fluid_cell.temperature;
0199   fluid_cell = bulk_info.get(tau_low, x_right, y_right, eta_left);
0200   cube[0][1][1][0] = fluid_cell.temperature;
0201   fluid_cell = bulk_info.get(tau_low, x_right, y_right, eta_right);
0202   cube[0][1][1][1] = fluid_cell.temperature;
0203   fluid_cell = bulk_info.get(tau_high, x_left, y_left, eta_left);
0204   cube[1][0][0][0] = fluid_cell.temperature;
0205   fluid_cell = bulk_info.get(tau_high, x_left, y_left, eta_right);
0206   cube[1][0][0][1] = fluid_cell.temperature;
0207   fluid_cell = bulk_info.get(tau_high, x_left, y_right, eta_left);
0208   cube[1][0][1][0] = fluid_cell.temperature;
0209   fluid_cell = bulk_info.get(tau_high, x_left, y_right, eta_right);
0210   cube[1][0][1][1] = fluid_cell.temperature;
0211   fluid_cell = bulk_info.get(tau_high, x_right, y_left, eta_left);
0212   cube[1][1][0][0] = fluid_cell.temperature;
0213   fluid_cell = bulk_info.get(tau_high, x_right, y_left, eta_right);
0214   cube[1][1][0][1] = fluid_cell.temperature;
0215   fluid_cell = bulk_info.get(tau_high, x_right, y_right, eta_left);
0216   cube[1][1][1][0] = fluid_cell.temperature;
0217   fluid_cell = bulk_info.get(tau_high, x_right, y_right, eta_right);
0218   cube[1][1][1][1] = fluid_cell.temperature;
0219
0220   if ((T_cut - cube[0][0][0][0]) * (cube[1][1][1][1] - T_cut) < 0.0)
0221     if ((T_cut - cube[0][0][1][1]) * (cube[1][1][0][0] - T_cut) < 0.0)
0222       if ((T_cut - cube[0][1][0][1]) * (cube[1][0][1][0] - T_cut) < 0.0)
0223         if ((T_cut - cube[0][1][1][0]) * (cube[1][0][0][1] - T_cut) < 0.0)
0224           if ((T_cut - cube[0][0][0][1]) * (cube[1][1][1][0] - T_cut) < 0.0)
0225             if ((T_cut - cube[0][0][1][0]) * (cube[1][1][0][1] - T_cut) < 0.0)
0226               if ((T_cut - cube[0][1][0][0]) * (cube[1][0][1][1] - T_cut) < 0.0)
0227                 if ((T_cut - cube[0][1][1][1]) * (cube[1][0][0][0] - T_cut) <
0228                     0.0)
0229                   intersect = false;
0230
0231   return (intersect);
0232 }
0233
0234 void SurfaceFinder::Find_full_hypersurface_4D() {
0235   auto grid_tau0 = bulk_info.Tau0();
0236   auto grid_tauf = bulk_info.TauMax();
0237   auto grid_x0 = bulk_info.XMin();
0238   auto grid_y0 = bulk_info.YMin();
0239   ;
0240   auto grid_eta0 = bulk_info.EtaMin();
0241   ;
0242
0243   Jetscape::real grid_dt = 0.1;
0244   Jetscape::real grid_dx = 0.2;
0245   Jetscape::real grid_dy = 0.2;
0246   Jetscape::real grid_deta = 0.2;
0247
0248   const int dim = 4;
0249   double lattice_spacing[dim];
0250   lattice_spacing[0] = grid_dt;
0251   lattice_spacing[1] = grid_dx;
0252   lattice_spacing[2] = grid_dy;
0253   lattice_spacing[3] = grid_deta;
0254
0255   std::unique_ptr<Cornelius> cornelius_ptr(new Cornelius());
0256   cornelius_ptr->init(dim, T_cut, lattice_spacing);
0257
0258   const int ntime = static_cast<int>((grid_tauf - grid_tau0) / grid_dt);
0259   const int nx = static_cast<int>(std::abs(2. * grid_x0) / grid_dx);
0260   const int ny = static_cast<int>(std::abs(2. * grid_y0) / grid_dy);
0261   const int neta = static_cast<int>(std::abs(2. * grid_eta0) / grid_deta);
0262
0263   double ****cube = new double ***[2];
0264   for (int i = 0; i < 2; i++) {
0265     cube[i] = new double **[2];
0266     for (int j = 0; j < 2; j++) {
0267       cube[i][j] = new double *[2];
0268       for (int k = 0; k < 2; k++) {
0269         cube[i][j][k] = new double[2];
0270         for (int l = 0; l < 2; l++) {
0271           cube[i][j][k][l] = 0.0;
0272         }
0273       }
0274     }
0275   }
0276
0277   for (int itime = 0; itime < ntime; itime++) {
0278     // loop over time evolution
0279     auto tau_local = grid_tau0 + (itime + 0.5) * grid_dt;
0280     for (int l = 0; l < neta; l++) {
0281       auto eta_local = grid_eta0 + (l + 0.5) * grid_deta;
0282       // loops over the transverse plane
0283       for (int i = 0; i < nx; i++) {
0284         // loops over the transverse plane
0285         auto x_local = grid_x0 + (i + 0.5) * grid_dx;
0286         for (int j = 0; j < ny; j++) {
0287           auto y_local = grid_y0 + (j + 0.5) * grid_dy;
0288           bool intersect =
0289               check_intersect_4D(tau_local, x_local, y_local, eta_local,
0290                                  grid_dt, grid_dx, grid_dy, grid_deta, cube);
0291           if (intersect) {
0292             cornelius_ptr->find_surface_4d(cube);
0293             for (int isurf = 0; isurf < cornelius_ptr->get_Nelements();
0294                  isurf++) {
0295               auto tau_center = (cornelius_ptr->get_centroid_elem(isurf, 0) +
0296                                  tau_local - grid_dt / 2.);
0297               auto x_center = (cornelius_ptr->get_centroid_elem(isurf, 1) +
0298                                x_local - grid_dx / 2.);
0299               auto y_center = (cornelius_ptr->get_centroid_elem(isurf, 2) +
0300                                y_local - grid_dy / 2.);
0301               auto eta_center = (cornelius_ptr->get_centroid_elem(isurf, 3) +
0302                                  eta_local - grid_deta / 2.);
0303
0304               auto da_tau = (cornelius_ptr->get_normal_elem(isurf, 0));
0305               auto da_x = (cornelius_ptr->get_normal_elem(isurf, 1));
0306               auto da_y = (cornelius_ptr->get_normal_elem(isurf, 2));
0307               auto da_eta = (cornelius_ptr->get_normal_elem(isurf, 3));
0308
0309               auto fluid_cell =
0310                   bulk_info.get(tau_center, x_center, y_center, eta_center);
0311               auto surface_cell = PrepareASurfaceCell(
0312                   tau_center, x_center, y_center, eta_center, da_tau, da_x,
0313                   da_y, da_eta, fluid_cell);
0314               surface_cell_list.push_back(surface_cell);
0315             }
0316           }
0317         }
0318       }
0319     }
0320   }
0321
0322   for (int i = 0; i < 2; i++) {
0323     for (int j = 0; j < 2; j++) {
0324       for (int k = 0; k < 2; k++) {
0325         delete[] cube[i][j][k];
0326       }
0327       delete[] cube[i][j];
0328     }
0329     delete[] cube[i];
0330   }
0331   delete[] cube;
0332 }
0333
0334 SurfaceCellInfo SurfaceFinder::PrepareASurfaceCell(
0335     Jetscape::real tau, Jetscape::real x, Jetscape::real y, Jetscape::real eta,
0336     Jetscape::real da0, Jetscape::real da1, Jetscape::real da2,
0337     Jetscape::real da3, const FluidCellInfo fluid_cell) {
0338
0339   SurfaceCellInfo temp_cell;
0340   temp_cell.tau = tau;
0341   temp_cell.x = x;
0342   temp_cell.y = y;
0343   temp_cell.eta = eta;
0344   temp_cell.d3sigma_mu[0] = da0;
0345   temp_cell.d3sigma_mu[1] = da1;
0346   temp_cell.d3sigma_mu[2] = da2;
0347   temp_cell.d3sigma_mu[3] = da3;
0348
0349   temp_cell.energy_density = fluid_cell.energy_density;
0350   temp_cell.entropy_density = fluid_cell.entropy_density;
0351   temp_cell.temperature = fluid_cell.temperature;
0352   temp_cell.pressure = fluid_cell.pressure;
0353   temp_cell.qgp_fraction = fluid_cell.qgp_fraction;
0354   temp_cell.mu_B = fluid_cell.mu_B;
0355   temp_cell.mu_C = fluid_cell.mu_C;
0356   temp_cell.mu_S = fluid_cell.mu_S;
0357
0358   temp_cell.vx = fluid_cell.vx;
0359   temp_cell.vy = fluid_cell.vy;
0360   temp_cell.vz = fluid_cell.vz;
0361
0362   for (int i = 0; i < 4; i++) {
0363     for (int j = 0; j < 4; j++) {
0364       temp_cell.pi[i][j] = fluid_cell.pi[i][j];
0365     }
0366   }
0367   temp_cell.bulk_Pi = fluid_cell.bulk_Pi;
0368
0369   return (temp_cell);
0370 }
0371
0372 } // namespace Jetscape