Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SurfaceFinder.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SurfaceFinder.cpp
1 // Copyright @ Chun Shen 2014
2 #include <iostream>
3 #include <sstream>
4 #include <fstream>
5 #include <cmath>
6 #include <iomanip>
7 #include <string>
8 #include <fstream>
9 
10 #include "./SurfaceFinder.h"
11 #include "./cornelius.h"
12 
13 using namespace std;
14 
15 SurfaceFinder::SurfaceFinder(void* hydroinfo_ptr_in,
16  ParameterReader* paraRdr_in) {
17  paraRdr = paraRdr_in;
18  hydro_type = paraRdr->getVal("hydro_type");
19  if (hydro_type == 0) {
20 #ifdef USE_HDF5
21  hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
22 #endif
23  } else {
24  hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
25  }
26  T_cut = paraRdr->getVal("T_cut");
27 }
28 
29 SurfaceFinder::SurfaceFinder(void* hydroinfo_ptr_in,
30  ParameterReader* paraRdr_in, double T_cut_in) {
31  paraRdr = paraRdr_in;
32  hydro_type = paraRdr->getVal("hydro_type");
33  if (hydro_type == 0) {
34 #ifdef USE_HDF5
35  hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
36 #endif
37  } else {
38  hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
39  }
40  T_cut = T_cut_in;
41 }
42 
43 SurfaceFinder::~SurfaceFinder() {}
44 
45 bool SurfaceFinder::check_intersect(double T_cut, double tau, double x,
46  double y, double dt, double dx, double dy,
47  double ***cube) {
48  hydrofluidCell *fluidCellptr = new hydrofluidCell();
49  bool intersect = true;
50 
51  double tau_low = tau - dt/2.;
52  double tau_high = tau + dt/2.;
53  double x_left = x - dx/2.;
54  double x_right = x + dx/2.;
55  double y_left = y - dy/2.;
56  double y_right = y + dy/2.;
57 
58  if (hydro_type == 0) {
59 #ifdef USE_HDF5
60  hydroinfo_ptr->getHydroinfo(tau_low, x_left, y_left, fluidCellptr);
61 #endif
62  } else {
63  hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_left, 0.0, tau_low,
64  fluidCellptr);
65  }
66  cube[0][0][0] = fluidCellptr->temperature;
67  if (hydro_type == 0) {
68 #ifdef USE_HDF5
69  hydroinfo_ptr->getHydroinfo(tau_low, x_left, y_right, fluidCellptr);
70 #endif
71  } else {
72  hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_right, 0.0, tau_low,
73  fluidCellptr);
74  }
75  cube[0][0][1] = fluidCellptr->temperature;
76  if (hydro_type == 0) {
77 #ifdef USE_HDF5
78  hydroinfo_ptr->getHydroinfo(tau_low, x_right, y_left, fluidCellptr);
79 #endif
80  } else {
81  hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_left, 0.0, tau_low,
82  fluidCellptr);
83  }
84  cube[0][1][0] = fluidCellptr->temperature;
85  if (hydro_type == 0) {
86 #ifdef USE_HDF5
87  hydroinfo_ptr->getHydroinfo(tau_low, x_right, y_right, fluidCellptr);
88 #endif
89  } else {
90  hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_right, 0.0, tau_low,
91  fluidCellptr);
92  }
93  cube[0][1][1] = fluidCellptr->temperature;
94  if (hydro_type == 0) {
95 #ifdef USE_HDF5
96  hydroinfo_ptr->getHydroinfo(tau_high, x_left, y_left, fluidCellptr);
97 #endif
98  } else {
99  hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_left, 0.0, tau_high,
100  fluidCellptr);
101  }
102  cube[1][0][0] = fluidCellptr->temperature;
103  if (hydro_type == 0) {
104 #ifdef USE_HDF5
105  hydroinfo_ptr->getHydroinfo(tau_high, x_left, y_right, fluidCellptr);
106 #endif
107  } else {
108  hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_right, 0.0, tau_high,
109  fluidCellptr);
110  }
111  cube[1][0][1] = fluidCellptr->temperature;
112  if (hydro_type == 0) {
113 #ifdef USE_HDF5
114  hydroinfo_ptr->getHydroinfo(tau_high, x_right, y_left, fluidCellptr);
115 #endif
116  } else {
117  hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_left, 0.0, tau_high,
118  fluidCellptr);
119  }
120  cube[1][1][0] = fluidCellptr->temperature;
121  if (hydro_type == 0) {
122 #ifdef USE_HDF5
123  hydroinfo_ptr->getHydroinfo(tau_high, x_right, y_right, fluidCellptr);
124 #endif
125  } else {
126  hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_right, 0.0, tau_high,
127  fluidCellptr);
128  }
129  cube[1][1][1] = fluidCellptr->temperature;
130 
131  if ((T_cut - cube[0][0][0])*(cube[1][1][1] - T_cut) < 0.0)
132  if ((T_cut - cube[0][1][0])*(cube[1][0][1] - T_cut) < 0.0)
133  if ((T_cut - cube[0][1][1])*(cube[1][0][0] - T_cut) < 0.0)
134  if ((T_cut - cube[0][0][1])*(cube[1][1][0] - T_cut) < 0.0)
135  intersect = false;
136 
137  delete fluidCellptr;
138  return(intersect);
139 }
140 
141 int SurfaceFinder::Find_full_hypersurface() {
142  ofstream output;
143  output.open("hyper_surface_2+1d.dat");
144 
145  double grid_tau0 = 0.0;
146  double grid_tauf = 0.0;
147  double grid_x0 = 0.0;
148  double grid_y0 = 0.0;
149  if (hydro_type == 1) {
150 #ifdef USE_HDF5
151  grid_tau0 = hydroinfo_ptr->getHydrogridTau0();
152  grid_tauf = hydroinfo_ptr->getHydrogridTaumax();
153  grid_x0 = hydroinfo_ptr->getHydrogridX0();
154  grid_y0 = hydroinfo_ptr->getHydrogridY0();
155 #endif
156  } else {
157  grid_tau0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
158  grid_tauf = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
159  grid_x0 = (- hydroinfo_MUSIC_ptr->get_hydro_x_max()
160  + hydroinfo_MUSIC_ptr->get_hydro_dx());
161  grid_y0 = grid_x0;
162  }
163 
164  double grid_dt = paraRdr->getVal("grid_dt");
165  double grid_dx = paraRdr->getVal("grid_dx");
166  double grid_dy = paraRdr->getVal("grid_dy");
167 
168  int dim = 3;
169  double *lattice_spacing = new double [dim];
170  lattice_spacing[0] = grid_dt;
171  lattice_spacing[1] = grid_dx;
172  lattice_spacing[2] = grid_dy;
173 
174  Cornelius* cornelius_ptr = new Cornelius();
175  cornelius_ptr->init(dim, T_cut, lattice_spacing);
176 
177  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt);
178  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx);
179  int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy);
180 
181  double ***cube = new double** [2];
182  for (int i = 0; i < 2; i++) {
183  cube[i] = new double* [2];
184  for (int j = 0; j < 2; j++) {
185  cube[i][j] = new double [2];
186  for (int k = 0; k < 2; k++)
187  cube[i][j][k] = 0.0;
188  }
189  }
190 
191  hydrofluidCell *fluidCellptr = new hydrofluidCell();
192 
193  for (int itime = 0; itime < ntime; itime++) {
194  // loop over time evolution
195  double tau_local = grid_tau0 + (itime + 0.5)*grid_dt;
196  for (int i = 0; i < nx; i++) {
197  // loops over the transverse plane
198  double x_local = grid_x0 + (i + 0.5)*grid_dx;
199  for (int j = 0; j < ny; j++) {
200  double y_local = grid_y0 + (j + 0.5)*grid_dy;
201  bool intersect = check_intersect(T_cut, tau_local, x_local,
202  y_local, grid_dt, grid_dx,
203  grid_dy, cube);
204  if (intersect) {
205  cornelius_ptr->find_surface_3d(cube);
206  for (int isurf = 0; isurf < cornelius_ptr->get_Nelements();
207  isurf++) {
208  double tau_center = (
209  cornelius_ptr->get_centroid_elem(isurf, 0)
210  + tau_local - grid_dt/2.);
211  double x_center = (
212  cornelius_ptr->get_centroid_elem(isurf, 1)
213  + x_local - grid_dx/2.);
214  double y_center = (
215  cornelius_ptr->get_centroid_elem(isurf, 2)
216  + y_local - grid_dy/2.);
217 
218  double da_tau =
219  cornelius_ptr->get_normal_elem(isurf, 0);
220  double da_x = cornelius_ptr->get_normal_elem(isurf, 1);
221  double da_y = cornelius_ptr->get_normal_elem(isurf, 2);
222 
223  if (hydro_type == 1) {
224 #ifdef USE_HDF5
225  hydroinfo_ptr->getHydroinfo(
226  tau_center, x_center, y_center, fluidCellptr);
227 #endif
228  } else {
229  hydroinfo_MUSIC_ptr->getHydroValues(
230  x_center, y_center, 0.0, tau_center,
231  fluidCellptr);
232  }
233 
234  output << scientific << setw(18) << setprecision(8)
235  << tau_center << " " << x_center << " "
236  << y_center << " "
237  << da_tau << " " << da_x << " "
238  << da_y << " "
239  << fluidCellptr->temperature << " "
240  << fluidCellptr->vx << " " << fluidCellptr->vy
241  << endl;
242 
243  }
244  }
245  }
246  }
247  }
248  output.close();
249 
250  delete fluidCellptr;
251  delete cornelius_ptr;
252  delete [] lattice_spacing;
253  for (int i = 0; i < 2; i++) {
254  for (int j = 0; j < 2; j++)
255  delete [] cube[i][j];
256  delete [] cube[i];
257  }
258  delete [] cube;
259  return 0;
260 }