6. Additional topics

6.1. Writing full spacetime evolution to disk

For certain applications, such as when calculating parton energy loss or photon production, it can be necessary to output the entire spacetime information of the hydrodynamic fields (\(\epsilon(X)\), \(u^\mu(X)\), …) to disk.

Writing this information to disk can create very large files, especially for 3+1D hydrodynamic simulations, and it should only be done when required.

Function “Grid_info::OutputEvolutionDataXYEta()” in source file “grid_info.cpp” provides an example of how to write the entire evolution of the plasma to disk.

void Grid_info::OutputEvolutionDataXYEta(...) {
    Util *util = new Util();

    int nx = DATA->nx;
    int ny = DATA->ny;
    int neta = DATA->neta;

    // MPI send and receive
    int sizeOfData = (nx+1)*(ny+1)*(neta);
    int to, from;

    [Code to take care of communication between parallel MPI instances]

    if (rank == 0) {
        [...]


        // set up file name
        const string out_name_xyeta = "evolution_xyeta.dat";
        const string out_name_W_xyeta = "evolution_Wmunu_over_epsilon_plus_P_xyeta.dat";
        const string out_name_bulk_xyeta = "evolution_bulk_pressure_xyeta.dat";
        string out_open_mode;
        FILE *out_file_xyeta;
        FILE *out_file_W_xyeta;
        FILE *out_file_bulk_xyeta;

        // If it's the first timestep, overwrite the previous file
        if (tau == DATA->tau0) {
           out_open_mode = "w";
        }
        else {
           out_open_mode = "a";
        }
        //If we output in binary, set the mode accordingly
        if (0 == DATA->outputBinaryEvolution) {
          out_open_mode += "b";
        }

        if (DATA->outputEvolutionData == 1) {
            out_file_xyeta =
                fopen(out_name_xyeta.c_str(), out_open_mode.c_str());
            out_file_W_xyeta =
                fopen(out_name_W_xyeta.c_str(), out_open_mode.c_str());
            out_file_bulk_xyeta =
                fopen(out_name_bulk_xyeta.c_str(), out_open_mode.c_str());
        }

        // Although it is a little confusing,
        // it is easiest to use (tau, x, y, eta) coordinates
        // and save at these points vx, vy, and vz. -CFY 11/16/2010
        float T1, u01, ux1, uy1, uz1, ueta1, utau1, epsilon1, rhob1, QGPfrac1;
        float entropy1;
        double eta;
        float Wtt, Wtx, Wty, Wtz, Wzz, Wxz, Wyz; // added by Maxime
        float Wtautau, Wtaux, Wtauy, Wtaueta, Wxx, Wxy, Wxeta, Wyy;
        float Wyeta, Wetaeta; // added by Maxime
        float div_factor, pressure1;

        int n_eta_step = DATA_ptr->output_evolution_every_N_eta;
        int n_x_step = DATA_ptr->output_evolution_every_N_x;
        int n_y_step = DATA_ptr->output_evolution_every_N_y;

        // No interpolation is necessary here!
        for (int ieta = 0; ieta < (DATA->neta)*size; ieta += n_eta_step) {
            // Recall that the eta grid is shifted by one cell
            // in the negative direction
            // For a boost-invariant system,
            // we want to output this information at eta=0,
            // not eta=-(DATA->eta_size)/2.0
            if (DATA->boost_invariant) {
                eta = 0.0;
            } else {
                eta = ((double)ieta)*(DATA->delta_eta)-(DATA->eta_size)/2.0;
            }

            for (int iy = 0; iy <= DATA->ny; iy += n_y_step) {
                for (int ix = 0; ix <= DATA->nx; ix += n_x_step) {
                        epsilon1 = epsFrom[ix][iy][ieta];
                        pressure1 = pFrom[ix][iy][ieta];
                        rhob1 = rhobFrom[ix][iy][ieta];
                        utau1 = utauFrom[ix][iy][ieta];
                        ux1 = uxFrom[ix][iy][ieta];
                        uy1 = uyFrom[ix][iy][ieta];
                        ueta1 = uetaFrom[ix][iy][ieta];

                        [...]

                    if (DATA->outputEvolutionData == 1) {
                        // exclude the actual coordinates from the output to save space:
                        if (0 == DATA->outputBinaryEvolution) {
                                fprintf(out_file_xyeta,"%e %e %e %e %e\n", T1*hbarc, QGPfrac1, ux1, uy1, uz1);
                                if (DATA->viscosity_flag) {
                                        if (1 == DATA->turn_on_shear) {
                                            fprintf(out_file_W_xyeta,"%e %e %e %e %e %e %e %e %e %e\n",Wtt,Wtx,Wty,Wtz,Wxx,Wxy,Wxz,Wyy,Wyz,Wzz);
                                        }
                                        if (1 == DATA->turn_on_bulk) {
                                            fprintf(out_file_bulk_xyeta,"%e %e %e\n", local_bulk, epsilon1+pressure1, cs2);
                                        }
                                }
                            } else {
                                float array[]={T1*hbarc, QGPfrac1, ux1, uy1, uz1};
                                fwrite(array,sizeof(float),5,out_file_xyeta);
                                // Write Wmunu/shear only if viscosity is on
                            // --- no need to fill a file with zeros in the ideal case
                                if (DATA->viscosity_flag) {
                                        if (1 == DATA->turn_on_shear) {
                                                float array2[]={Wtt,Wtx,Wty,Wtz,Wxx,Wxy,Wxz,Wyy,Wyz,Wzz};
                                                fwrite(array2,sizeof(float),10,out_file_W_xyeta);
                                        }
                                        if (1 == DATA->turn_on_bulk) {
                                                float array2[]={local_bulk, epsilon1+pressure1, cs2};
                                            fwrite(array2,sizeof(float),3,out_file_bulk_xyeta);
                                        }
                                }
                            }
                    }
                }/* ix */
            }/* iy */
        }/* ieta */
        if (DATA->outputEvolutionData == 1) {
            fclose(out_file_xyeta);
            fclose(out_file_W_xyeta);
            fclose(out_file_bulk_xyeta);
        }

        [...]
    }
    [...]

}/* OutputEvolutionDataXYEta */