external_packages/hydro_from_external_file/SurfaceFinder.cpp

0001 // Copyright @ Chun Shen 2014
0002 #include <iostream>
0003 #include <sstream>
0004 #include <fstream>
0005 #include <cmath>
0006 #include <iomanip>
0007 #include <string>
0008 #include <fstream>
0009
0010 #include "./SurfaceFinder.h"
0011 #include "./cornelius.h"
0012
0013 using namespace std;
0014
0015 SurfaceFinder::SurfaceFinder(void* hydroinfo_ptr_in,
0016                              ParameterReader* paraRdr_in) {
0017     paraRdr = paraRdr_in;
0018     hydro_type = paraRdr->getVal("hydro_type");
0019     if (hydro_type == 0) {
0020 #ifdef USE_HDF5
0021         hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
0022 #endif
0023     } else {
0024         hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
0025     }
0026     T_cut = paraRdr->getVal("T_cut");
0027 }
0028
0029 SurfaceFinder::SurfaceFinder(void* hydroinfo_ptr_in,
0030                              ParameterReader* paraRdr_in, double T_cut_in) {
0031     paraRdr = paraRdr_in;
0032     hydro_type = paraRdr->getVal("hydro_type");
0033     if (hydro_type == 0) {
0034 #ifdef USE_HDF5
0035         hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
0036 #endif
0037     } else {
0038         hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
0039     }
0040     T_cut = T_cut_in;
0041 }
0042
0043 SurfaceFinder::~SurfaceFinder() {}
0044
0045 bool SurfaceFinder::check_intersect(double T_cut, double tau, double x,
0046                                     double y, double dt, double dx, double dy,
0047                                     double ***cube) {
0048     hydrofluidCell *fluidCellptr = new hydrofluidCell();
0049     bool intersect = true;
0050
0051     double tau_low = tau - dt/2.;
0052     double tau_high = tau + dt/2.;
0053     double x_left = x - dx/2.;
0054     double x_right = x + dx/2.;
0055     double y_left = y - dy/2.;
0056     double y_right = y + dy/2.;
0057
0058     if (hydro_type == 0) {
0059 #ifdef USE_HDF5
0060         hydroinfo_ptr->getHydroinfo(tau_low, x_left, y_left, fluidCellptr);
0061 #endif
0062     } else {
0063         hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_left, 0.0, tau_low,
0064                                             fluidCellptr);
0065     }
0066     cube[0][0][0] = fluidCellptr->temperature;
0067     if (hydro_type == 0) {
0068 #ifdef USE_HDF5
0069         hydroinfo_ptr->getHydroinfo(tau_low, x_left, y_right, fluidCellptr);
0070 #endif
0071     } else {
0072         hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_right, 0.0, tau_low,
0073                                             fluidCellptr);
0074     }
0075     cube[0][0][1] = fluidCellptr->temperature;
0076     if (hydro_type == 0) {
0077 #ifdef USE_HDF5
0078         hydroinfo_ptr->getHydroinfo(tau_low, x_right, y_left, fluidCellptr);
0079 #endif
0080     } else {
0081         hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_left, 0.0, tau_low,
0082                                             fluidCellptr);
0083     }
0084     cube[0][1][0] = fluidCellptr->temperature;
0085     if (hydro_type == 0) {
0086 #ifdef USE_HDF5
0087         hydroinfo_ptr->getHydroinfo(tau_low, x_right, y_right, fluidCellptr);
0088 #endif
0089     } else {
0090         hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_right, 0.0, tau_low,
0091                                             fluidCellptr);
0092     }
0093     cube[0][1][1] = fluidCellptr->temperature;
0094     if (hydro_type == 0) {
0095 #ifdef USE_HDF5
0096         hydroinfo_ptr->getHydroinfo(tau_high, x_left, y_left, fluidCellptr);
0097 #endif
0098     } else {
0099         hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_left, 0.0, tau_high,
0100                                             fluidCellptr);
0101     }
0102     cube[1][0][0] = fluidCellptr->temperature;
0103     if (hydro_type == 0) {
0104 #ifdef USE_HDF5
0105         hydroinfo_ptr->getHydroinfo(tau_high, x_left, y_right, fluidCellptr);
0106 #endif
0107     } else {
0108         hydroinfo_MUSIC_ptr->getHydroValues(x_left, y_right, 0.0, tau_high,
0109                                             fluidCellptr);
0110     }
0111     cube[1][0][1] = fluidCellptr->temperature;
0112     if (hydro_type == 0) {
0113 #ifdef USE_HDF5
0114         hydroinfo_ptr->getHydroinfo(tau_high, x_right, y_left, fluidCellptr);
0115 #endif
0116     } else {
0117         hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_left, 0.0, tau_high,
0118                                             fluidCellptr);
0119     }
0120     cube[1][1][0] = fluidCellptr->temperature;
0121     if (hydro_type == 0) {
0122 #ifdef USE_HDF5
0123         hydroinfo_ptr->getHydroinfo(tau_high, x_right, y_right, fluidCellptr);
0124 #endif
0125     } else {
0126         hydroinfo_MUSIC_ptr->getHydroValues(x_right, y_right, 0.0, tau_high,
0127                                             fluidCellptr);
0128     }
0129     cube[1][1][1] = fluidCellptr->temperature;
0130
0131     if ((T_cut - cube[0][0][0])*(cube[1][1][1] - T_cut) < 0.0)
0132         if ((T_cut - cube[0][1][0])*(cube[1][0][1] - T_cut) < 0.0)
0133             if ((T_cut - cube[0][1][1])*(cube[1][0][0] - T_cut) < 0.0)
0134                 if ((T_cut - cube[0][0][1])*(cube[1][1][0] - T_cut) < 0.0)
0135                     intersect = false;
0136
0137     delete fluidCellptr;
0138     return(intersect);
0139 }
0140
0141 int SurfaceFinder::Find_full_hypersurface() {
0142     ofstream output;
0143     output.open("hyper_surface_2+1d.dat");
0144
0145     double grid_tau0 = 0.0;
0146     double grid_tauf = 0.0;
0147     double grid_x0 = 0.0;
0148     double grid_y0 = 0.0;
0149     if (hydro_type == 1) {
0150 #ifdef USE_HDF5
0151         grid_tau0 = hydroinfo_ptr->getHydrogridTau0();
0152         grid_tauf = hydroinfo_ptr->getHydrogridTaumax();
0153         grid_x0 = hydroinfo_ptr->getHydrogridX0();
0154         grid_y0 = hydroinfo_ptr->getHydrogridY0();
0155 #endif
0156     } else {
0157         grid_tau0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
0158         grid_tauf = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
0159         grid_x0 = (- hydroinfo_MUSIC_ptr->get_hydro_x_max()
0160                    + hydroinfo_MUSIC_ptr->get_hydro_dx());
0161         grid_y0 = grid_x0;
0162     }
0163
0164     double grid_dt = paraRdr->getVal("grid_dt");
0165     double grid_dx = paraRdr->getVal("grid_dx");
0166     double grid_dy = paraRdr->getVal("grid_dy");
0167
0168     int dim = 3;
0169     double *lattice_spacing = new double [dim];
0170     lattice_spacing[0] = grid_dt;
0171     lattice_spacing[1] = grid_dx;
0172     lattice_spacing[2] = grid_dy;
0173
0174     Cornelius* cornelius_ptr = new Cornelius();
0175     cornelius_ptr->init(dim, T_cut, lattice_spacing);
0176
0177     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt);
0178     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx);
0179     int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy);
0180
0181     double ***cube = new double** [2];
0182     for (int i = 0; i < 2; i++) {
0183         cube[i] = new double* [2];
0184         for (int j = 0; j < 2; j++) {
0185             cube[i][j] = new double [2];
0186             for (int k = 0; k < 2; k++)
0187                 cube[i][j][k] = 0.0;
0188         }
0189     }
0190
0191     hydrofluidCell *fluidCellptr = new hydrofluidCell();
0192
0193     for (int itime = 0; itime < ntime; itime++) {
0194         // loop over time evolution
0195         double tau_local = grid_tau0 + (itime + 0.5)*grid_dt;
0196         for (int i = 0; i < nx; i++) {
0197             // loops over the transverse plane
0198             double x_local = grid_x0 + (i + 0.5)*grid_dx;
0199             for (int j = 0; j < ny; j++) {
0200                 double y_local = grid_y0 + (j + 0.5)*grid_dy;
0201                 bool intersect = check_intersect(T_cut, tau_local, x_local,
0202                                                  y_local, grid_dt, grid_dx,
0203                                                  grid_dy, cube);
0204                 if (intersect) {
0205                     cornelius_ptr->find_surface_3d(cube);
0206                     for (int isurf = 0; isurf < cornelius_ptr->get_Nelements();
0207                          isurf++) {
0208                         double tau_center = (
0209                                 cornelius_ptr->get_centroid_elem(isurf, 0)
0210                                 + tau_local - grid_dt/2.);
0211                         double x_center = (
0212                                 cornelius_ptr->get_centroid_elem(isurf, 1)
0213                                 + x_local - grid_dx/2.);
0214                         double y_center = (
0215                                 cornelius_ptr->get_centroid_elem(isurf, 2)
0216                                 + y_local - grid_dy/2.);
0217
0218                         double da_tau =
0219                                 cornelius_ptr->get_normal_elem(isurf, 0);
0220                         double da_x = cornelius_ptr->get_normal_elem(isurf, 1);
0221                         double da_y = cornelius_ptr->get_normal_elem(isurf, 2);
0222
0223                         if (hydro_type == 1) {
0224 #ifdef USE_HDF5
0225                             hydroinfo_ptr->getHydroinfo(
0226                                 tau_center, x_center, y_center, fluidCellptr);
0227 #endif
0228                         } else {
0229                             hydroinfo_MUSIC_ptr->getHydroValues(
0230                                 x_center, y_center, 0.0, tau_center,
0231                                 fluidCellptr);
0232                         }
0233
0234                         output << scientific << setw(18) << setprecision(8)
0235                                << tau_center << "   " << x_center << "   "
0236                                << y_center << "   "
0237                                << da_tau << "   " << da_x << "   "
0238                                << da_y << "   "
0239                                << fluidCellptr->temperature << "   "
0240                                << fluidCellptr->vx << "   " << fluidCellptr->vy
0241                                << endl;
0242
0243                     }
0244                 }
0245             }
0246         }
0247     }
0248     output.close();
0249
0250     delete fluidCellptr;
0251     delete cornelius_ptr;
0252     delete [] lattice_spacing;
0253     for (int i = 0; i < 2; i++) {
0254         for (int j = 0; j < 2; j++)
0255             delete [] cube[i][j];
0256         delete [] cube[i];
0257     }
0258     delete [] cube;
0259     return 0;
0260 }