/**
 * @file output.cpp
 * @author Christian Laguerre <christian.laguerre@math.cnrs.fr> (2012-2015)
 * @version 1.07.01
 * @date 2017-02-28
 *
 * @brief %Output format
 * @details 
 * Common part for all the output formats.
 *
 * @copyright License Cecill-V2 \n
 * <http://www.cecill.info/licences/Licence_CeCILL_V2-en.html>
 *
 * (c) CNRS - Universite d'Orleans - BRGM (France)
 */
/* 
 *
 * This file is part of FullSWOF_2D software. 
 * <https://sourcesup.renater.fr/projects/fullswof-2d/> 
 *
 * FullSWOF_2D = Full Shallow-Water equations for Overland Flow, 
 * in two dimensions of space.
 * This software is a computer program whose purpose is to compute
 * solutions for 2D Shallow-Water equations.
 *
 * LICENSE
 *
 * This software is governed by the CeCILL license under French law and
 * abiding by the rules of distribution of free software. You can use,
 * modify and/ or redistribute the software under the terms of the CeCILL
 * license as circulated by CEA, CNRS and INRIA at the following URL
 * <http://www.cecill.info>.
 *
 * As a counterpart to the access to the source code and rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty and the software's author, the holder of the
 * economic rights, and the successive licensors have only limited
 * liability.
 *
 * In this respect, the user's attention is drawn to the risks associated
 * with loading, using, modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean that it is complicated to manipulate, and that also
 * therefore means that it is reserved for developers and experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and, more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL license and that you accept its terms.
 *
 ******************************************************************************/


#include "output.hpp"

Output::Output(Parameters & par):NXCELL(par.get_Nxcell()),NYCELL(par.get_Nycell()),DX(par.get_dx()),DY(par.get_dy()){

  /**
   * @details  
   * Defines the names of the outputs.\n
   * If run in DEBUG mode, writes the header of the file 'boundaries_flux.dat', 'check_vol.dat', 'flux_boundaries_BT.dat' and 'flux_boundaries_LR.dat'.
   * @param[in] par parameter, contains all the values from the parameters file.
   * @warning Impossible to open the *** file. Verify if the directory *** exists.
   * @note If 'boundaries_flux.dat', 'check_vol.dat', 'flux_boundaries_BT.dat' or 'flux_boundaries_LR.dat' cannot be opened, the code will exit with failure termination code.
   */
  
  
  outputDirectory = par.get_outputDirectory();
  
  namefile_res = outputDirectory + "results.dat";
  
  namefile_init = outputDirectory + "huz_initial.dat";
  
  namefile_final = outputDirectory + "huz_final.dat";
  
#ifdef DEBUG
  
  namefile_check_volume = outputDirectory + "check_vol.dat";
  
  namefile_Bound_flux = outputDirectory + "boundaries_flux.dat";
  
  namefile_Bound_flux_BT= outputDirectory + "flux_boundaries_BT.dat";
  namefile_Bound_flux_LR= outputDirectory + "flux_boundaries_LR.dat";
  
  ofstream Boundaries_flux_f(namefile_Bound_flux.c_str(),ios::out);
  if (!Boundaries_flux_f){ 
    cerr << "impossible to open the "<< namefile_Bound_flux  << " file\n"; // if the name of the file is different, then impossible to open
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }

  Boundaries_flux_f.precision(10);

  Boundaries_flux_f << "################################################################## "<< endl;
  Boundaries_flux_f << "# Generated by "<< VERSION << endl; 
  Boundaries_flux_f << "################################################################## "<< endl;
  Boundaries_flux_f << "# Cumulated fluxes on the left, right, bottom and top boundaries at each time"<<endl;
  Boundaries_flux_f << "# time" << "\t"<< setw(9) << "left cumulative flux (m^2)" << "\t"<< setw(9) << "right cumulative flux (m^2)"  << "\t"<< setw(9) << "bottom cumulative flux (m^2)" << "\t"<<setw(9) << "top cumulative flux (m^2)" << endl;
  
  Boundaries_flux_f.close();
  
  
  ofstream check_vol_f(namefile_check_volume.c_str(),ios::out);
  if (!check_vol_f){ 
    cerr << "impossible to open the "<< namefile_check_volume << " file\n"; // if the name of the file is different, then impossible to open
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }

  check_vol_f.precision(10);

  check_vol_f << "################################################################## "<< endl;
  check_vol_f << "# Generated by "<< VERSION << endl; 
  check_vol_f << "################################################################## "<< endl;
  check_vol_f << "# time" << "\t"<< setw(9)  <<"Vol_of_tot" << "\t" <<setw(9) << "Vol_inf_tot" << "\t" <<setw(9) << "Vol_rain_tot"<< "\t" <<setw(9)<< "Vol_bound_tot"   << endl;
  
  check_vol_f.close();
  
  
  ofstream Boundaries_flux_BT_f(namefile_Bound_flux_BT.c_str(),ios::out);
  if (!Boundaries_flux_BT_f){
    cerr << "Impossible to open the " << namefile_Bound_flux_BT.c_str() <<" file\n";
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }

  Boundaries_flux_BT_f.precision(10);

  Boundaries_flux_BT_f << "################################################################## "<< endl;
  Boundaries_flux_BT_f << "# Generated by "<< VERSION << endl; 
  Boundaries_flux_BT_f << "################################################################## "<< endl;
  Boundaries_flux_BT_f << "# Flux on the bottom and top boundaries at each time"<<endl;
  Boundaries_flux_BT_f << "# (i-0.5)*dx" << "\t"<< setw(9) << "bottom flux (m^2/s)" << "\t"<< setw(9) << "top flux (m^2/s)" <<  "\t" << setw(9) << " time" << endl;
  Boundaries_flux_BT_f.close();
  
  ofstream Boundaries_flux_LR_f(namefile_Bound_flux_LR.c_str(),ios::out);
  if (!Boundaries_flux_LR_f){
    cerr << "Impossible to open the " << namefile_Bound_flux_LR.c_str() <<" file\n";
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }

  Boundaries_flux_LR_f.precision(10);

  Boundaries_flux_LR_f << "################################################################## "<< endl;
  Boundaries_flux_LR_f << "# Generated by "<< VERSION << endl;
  Boundaries_flux_LR_f << "################################################################## "<< endl;
  Boundaries_flux_LR_f << "# Flux on the left and right boundaries at each time"<<endl;  
  Boundaries_flux_LR_f << "# (j-0.5)*dy" << "\t"<< setw(9) << "left flux (m^2/s)" << "\t"<< setw(9) << "right flux (m^2/s)" <<  "\t" << setw(9) << " time" << endl;
  Boundaries_flux_LR_f .close();
#endif
  
  
  //These variables are necessary to compute the cumulated flux during the time
  fluxy0_cum_T = 0. ;
  fluxNycell_cum_T= 0. ;
  fluxx0_cum_T = 0. ;
  fluxNxcell_cum_T = 0. ;
  
}//end constructor


void Output:: initial(const  TAB & z, const TAB & h, const  TAB & u, const  TAB & v) const{
  
  /**
   * @details
   * @param[in] z topography.
   * @param[in] h water height.
   * @param[in] u first component of the velocity.
   * @param[in] v second component of the velocity.
   * @warning Impossible to open the *** file. Verify if the directory *** exists.
   * @note If huz_initial.dat cannot be opened, the code will exit with failure termination code.
   */
  
  ofstream initial(namefile_init.c_str(),ios::out);
  if (!initial){ 
    cerr << "impossible to open the "<< namefile_init << " file\n"; // if the name of the file is different, then impossible to open
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }
  initial << "################################################################## "<< endl;
  initial << "# Generated by "<< VERSION << endl; 
  initial << "################################################################## "<< endl;
  initial << "# (i-0.5)*dx" << "\t" << "(j-0.5)*dy"  << "\t" <<setw(9) << "h[i][j]" << "\t" <<setw(9) << "u[i][j]" << "\t" <<setw(9) << "v[i][j]"<< "\t" << setw(9) << "h[i][j]+z[i][j]" << "\t" << setw(9) << "z[i][j]" << endl;
  
  for (int i=1 ; i<=NXCELL ; i++){
    for (int j=1 ; j<=NYCELL ; j++){
      initial << (i-0.5)*DX << " \t" << (j-0.5)*DY  << "\t" <<setw(9) << h[i][j] << "\t" <<setw(9) <<u[i][j] << "\t" <<setw(9) << v[i][j]  << "\t" << setw(9) << h[i][j]+z[i][j] << "\t" << setw(9) << z[i][j] << endl;
      }// end of j loop
      initial << endl;
      
    }// end of i loop
    
    initial.close();
}


void Output:: check_vol(const SCALAR & time, const SCALAR & dt, const SCALAR & Vol_rain_tot, const SCALAR & Vol_inf, const SCALAR & Vol_of, const SCALAR & Vol_bound_tot) const {
  
  /**
   * @details
   * @param[in] time current time.
   * @param[in] dt time step (unused).
   * @param[in] Vol_rain_tot total rain volume.
   * @param[in] Vol_inf volume of infiltrated water.
   * @param[in] Vol_of volume of overland flow.
   * @param[in] Vol_bound_tot total volume of water at the boundary.
   */
  
  (void) dt; //unused variable
  
  ofstream check_vol_f(namefile_check_volume.c_str(),ios::app);
  check_vol_f << time << "\t"<< setw(9) << Vol_of  << "\t" <<setw(9) << Vol_inf << "\t" <<setw(9) << Vol_rain_tot<< "\t" <<setw(9)<< Vol_bound_tot << endl;
  
  check_vol_f.close();
  
}

void Output:: result(const SCALAR & time, const clock_t & cpu, const SCALAR & Vol_rain, const SCALAR & Vol_inf, const SCALAR & Vol_of , const SCALAR & FROUDE, const int & NBITER, const SCALAR & vol_output) const{
  
  /**
   * @details 
   * @param[in] time elapsed time.
   * @param[in] cpu CPU time.
   * @param[in] Vol_rain total rain volume.
   * @param[in] Vol_inf total volume of infiltrated water.
   * @param[in] Vol_of total volume of overland flow.
   * @param[in] FROUDE mean Froude number (in space) at the final time.
   * @param[in] NBITER number of time steps.
   * @param[in] vol_output total outflow volume at the boundary.
   * @warning Impossible to open the *** file. Verify if the directory *** exists.
   * @note If results.dat cannot be opened, the code will exit with failure termination code.
   */
  
  ofstream  result_f(namefile_res.c_str(),ios::out);
  if (!result_f){ 
    cerr << "impossible to open the "<< namefile_res << "\n"; // if the name of the file is different, then impossible to open
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }
  result_f << "################################################################## "<< endl;
  result_f << "# Generated by "<< VERSION << endl; 
  result_f << "################################################################## "<< endl;
  result_f << "Infiltrated volume: " << Vol_inf << endl;
  result_f << endl;
  result_f << "Stream volume: " << Vol_of << endl;
  result_f << endl;
  result_f << "Complete volume (Inf+Stream): " << Vol_of+Vol_inf << endl;
  result_f << endl;
  result_f << "Volume of the rain: " << Vol_rain << endl;
  result_f << endl;
  result_f << "Outflow volume at the boundaries: " << vol_output << endl;
  result_f << endl;
  result_f << "Time of the computation (seconds): " << time << endl;
  result_f << endl;
  result_f << "CPU time (clock ticks): "<< cpu<<endl;
  result_f << endl;
  
  //to give the number of iterations in the algorithm
  result_f << "Number of iterations in the algorithm: " << NBITER << endl;
  result_f << endl;
  
  
  result_f << "Froude number: " << FROUDE << endl;
  result_f << endl;
  
  result_f.close();
}


SCALAR Output::boundaries_flux(const SCALAR & time, const TAB & flux_u, const TAB & flux_v, const SCALAR & dt, const SCALAR & dt_first, const int & ORDER, const int & verif) {
  
  /**
   * @details
   * @param[in] time current time.
   * @param[in] flux_u flux on the left and right boundaries (m^2/s).
   * @param[in] flux_v flux on the bottom and top boundaries (m^2/s).
   * @param[in] dt current time step.
   * @param[in] dt_first previous time step.
   * @param[in] ORDER order of scheme.
   * @param[in] verif parameter to know if we removed the computation with the previous time step (dt_first).
   */
  
  SCALAR fluxNycell=0.;
  SCALAR fluxNxcell=0.;
  SCALAR fluxy0=0.;
  SCALAR fluxx0=0.;
  if (1 == verif){  
    for (int i=1 ; i<NXCELL+1 ; i++){
      fluxy0 +=  flux_v[i][1]; 
      fluxNycell += flux_v[i][NYCELL+1];
    }// end of i loop
    for (int j=1 ; j<NYCELL+1 ; j++){
      fluxx0 += flux_u[1][j];
      fluxNxcell += flux_u[NXCELL+1][j];
    }// end of j loop
  }
  
  fluxy0_cum_T =  fluxy0_cum_T+(fluxy0*(dt-dt_first*(1-verif))*(1./ORDER));
  fluxNycell_cum_T=  fluxNycell_cum_T+(fluxNycell*(dt-dt_first*(1-verif))*(1./ORDER));
  fluxx0_cum_T = fluxx0_cum_T+(fluxx0*(dt-dt_first*(1-verif))*(1./ORDER));
  fluxNxcell_cum_T = fluxNxcell_cum_T+(fluxNxcell*(dt-dt_first*(1-verif))*(1./ORDER));
  
#ifdef DEBUG
  
  ofstream Boundaries_flux_f(namefile_Bound_flux.c_str(),ios::app);
  Boundaries_flux_f << time << "\t"<< setw(9) << fluxx0_cum_T  << "\t"<< setw(9)<< fluxNxcell_cum_T  <<"\t"<< setw(9) << fluxy0_cum_T << "\t"<< setw(9)<< fluxNycell_cum_T  << endl;
  Boundaries_flux_f.close();
#endif
  
  return (DY*(fluxNxcell_cum_T-fluxx0_cum_T)+(fluxNycell_cum_T-fluxy0_cum_T)*DX);
}

void Output::boundaries_flux_LR(const SCALAR & time, const TAB & LR_flux) const{
  
  /**
   * @details
   * @param[in] time current time.
   * @param[in] LR_flux flux on the left and right boundaries (m^2/s).
   */
  
  ofstream Boundaries_flux_LR_f(namefile_Bound_flux_LR.c_str(),ios::app);
  Boundaries_flux_LR_f << "# time: "<<  time << " seconds "<<endl;
  for (int j=1 ; j<NYCELL+1 ; j++){
    Boundaries_flux_LR_f << (j-0.5)*DY << "\t"<< setw(9) << LR_flux[1][j] << "\t"<< setw(9)<<LR_flux[NXCELL+1][j] << "\t" << setw(9) << time << endl;
  }// end of j loop
  Boundaries_flux_LR_f << endl;
  Boundaries_flux_LR_f.close();
  
}

void Output::boundaries_flux_BT(const SCALAR & time, const TAB & BT_flux) const{
  
  /**
   * @details
   * @param[in] time current time.
   * @param[in] BT_flux flux on the bottom and tom boundaries (m^2/s).
   */
  
  ofstream Boundaries_flux_BT_f(namefile_Bound_flux_BT.c_str(),ios::app);
  Boundaries_flux_BT_f << "# time: "<<  time << " seconds "<<endl;
  for (int i=1 ; i<NXCELL+1 ; i++){
    Boundaries_flux_BT_f  << (i-0.5)*DX << "\t"<< setw(9) << BT_flux[i][1] << "\t"<< setw(9)<<BT_flux[i][NYCELL+1] << "\t" << setw(9) << time << endl;
  }// end of i loop
  Boundaries_flux_BT_f << endl;
  Boundaries_flux_BT_f.close();
  
}

void Output::final(const TAB & z, const TAB & h, const TAB & u, const TAB & v) const{
  
  /**
   * @details
   * If the water height is too small, we replace it by 0, the velocities and discharge are null and the Froude number does not exist.
   * @param[in] z topography.
   * @param[in] h water height.
   * @param[in] u first component of the velocity.
   * @param[in] v second component of the velocity.
   * @warning Impossible to open the *** file. Verify if the directory *** exists.
   * @note If huz_final.dat cannot be opened, the code will exit with failure termination code.
   */
  
  ofstream final(namefile_final.c_str(),ios::out);
  if (!final){ 
    cerr << "impossible to open the "<< namefile_final << " file\n"; // if the name of the file is different, then impossible to open
    cerr << "Verify if the directory " << outputDirectory <<" exists\n";
    exit(EXIT_FAILURE);
  }
  final << "################################################################## "<< endl;
  final << "# Generated by "<< VERSION << endl; 
  final << "################################################################## "<< endl;
  final << "# Final Time" << endl;
  final << "#(i-0.5)*dx" << "\t" << "(j-0.5)*dy" << "\t" << "h[i][j]" << "\t" << "u[i][j]" << "\t" << "v[i][j]"  << "\t" << "h[i][j]+z[i][j]" << "\t" << "z[i][j]" << "\t" << "velocity_norm[i][j]" << "\t" << "Fr[i][j] (Froude)" << "\t" << "qx[i][j]" << "\t" << "qy[i][j]" << "\t" << "q[i][j]" << endl;
  final << "#" << endl;
  for (int i=1 ; i<NXCELL+1 ; i++){
    for (int j=1 ; j<NYCELL+1 ; j++){
      if (h[i][j]<HE_CA){
        final << (i-0.5)*DX << "\t"<< setw(9) << (j-0.5)*DY << "\t"<< setw(9) << "0.0" << "\t"<< setw(9) << "0.0" << "\t"<< setw(9) << "0.0" << "\t"<< setw(9) << z[i][j] << "\t"<< setw(9) << z[i][j] << "\t"<< setw(9) << "0.0" << "\t"<< setw(9) << "NaN" << "\t"<< setw(9) << "0.0"  << "\t"<< setw(9) << "0.0" << "\t"<< setw(9) << "0.0" << endl;
      }
      else{
        final << (i-0.5)*DX << "\t"<< setw(9) << (j-0.5)*DY << "\t"<< setw(9) << h[i][j] << "\t"<< setw(9) << u[i][j] << "\t"<< setw(9) << v[i][j] << "\t"<< setw(9) << h[i][j]+z[i][j] << "\t"<< setw(9) << z[i][j] << "\t"<< setw(9) << sqrt(pow(u[i][j],2)+pow(v[i][j],2)) << "\t"<< setw(9) << max(sqrt(pow(u[i][j],2)+pow(v[i][j],2))/sqrt(GRAV*h[i][j]),0.) << "\t"<< setw(9) << h[i][j]*u[i][j]*DY  << "\t"<< setw(9) << h[i][j]*v[i][j]*DX << "\t"<< setw(9) << h[i][j]*sqrt(pow(u[i][j],2)+pow(v[i][j],2)) << endl;
      }
    }// end of j loop
    final << endl;
    
  }// end of i loop
  
  final.close();
}

Output::~Output(){
}