.. |Tmunu| replace:: :math:`T^{\mu\nu}`
.. |e| replace:: :math:`\epsilon`
.. |u| replace:: :math:`u^\mu`
.. |pimunu| replace:: :math:`\pi^{\mu\nu}`
.. |Pi| replace:: :math:`\Pi`
.. |Tmunu(X)| replace:: :math:`T^{\mu\nu}(X)`
.. |e(X)| replace:: :math:`\epsilon(X)`
.. |u(X)| replace:: :math:`u^\mu(X)`
.. |pimunu(X)| replace:: :math:`\pi^{\mu\nu}(X)`
.. |Pi(X)| replace:: :math:`\Pi(X)`

.. role:: bash(code)
   :language: bash

Additional topics
=====================

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 (|e(X)|, |u(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.

.. code-block :: c++

    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 */