/**
 * @file order2.cpp
 * @author Olivier Delestre <olivierdelestre41@yahoo.fr> (2008)
 * @author Christian Laguerre <christian.laguerre@math.cnrs.fr> (2012-2015)
 * @version 1.07.01
 * @date 2017-01-05
 *
 * @brief Order 2 scheme
 * @details 
 * Numerical scheme: at order 2.
 *
 * @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 "order2.hpp"

Order2::Order2(Parameters & par):Scheme(par){
  
  /**
   * @details   
   * Initializations, definition of the reconstruction and creation of 3 vectors for this reconstruction.
   * @param[in] par parameter, contains all the values from the parameters file.
   */
  
  Vin1.resize(NXCELL+1); // i : 1->NXCELL
  Vin2.resize(NXCELL+1); // i : 1->NXCELL
  
  hsa.resize(NXCELL+1); // i : 1->NXCELL
  usa.resize(NXCELL+1); // i : 1->NXCELL
  vsa.resize(NXCELL+1); // i : 1->NXCELL
  qsa1.resize(NXCELL+1); // i : 1->NXCELL
  qsa2.resize(NXCELL+1); // i : 1->NXCELL
  
  for (int i=1 ; i<=NXCELL ; i++){
    Vin1[i].resize(NYCELL+1); // l : 1->NYCELL
    Vin2[i].resize(NYCELL+1); // l : 1->NYCELL
    hsa[i].resize(NYCELL+1); // l : 1->NYCELL
    usa[i].resize(NYCELL+1); // l : 1->NYCELL
    vsa[i].resize(NYCELL+1); // l : 1->NYCELL
    qsa1[i].resize(NYCELL+1); // l : 1->NYCELL
    qsa2[i].resize(NYCELL+1); // l : 1->NYCELL
   
  }
  
  rec = new Choice_reconstruction(par,z);
  
  for (int i=1 ; i<=NXCELL ; i++){
    for (int j=1 ; j<=NYCELL ; j++){
      Vin1[i][j] = 0.;
      Vin2[i][j] = 0.;
    }
  }
}

Order2::~Order2(){
  
  for (int i=1 ; i<=NXCELL ; i++){
    Vin1[i].clear();
    Vin2[i].clear();
    hsa[i].clear();
    usa[i].clear();
    vsa[i].clear();
    qsa1[i].clear();
    qsa2[i].clear();
  }
  
  hsa.clear();
  usa.clear();
  vsa.clear();
  qsa1.clear();
  qsa2.clear();
  
  Vin1.clear();
  Vin2.clear();
  
  if (rec != NULL){
    delete rec;
    rec = NULL;
  }
}

void Order2::calcul(){
  
  /**
   * @details Performs the second order numerical scheme.
   * @note In DEBUG mode, the programme will save another file with volumes of water.
   * @warning order2: WARNING: the computation finished because the maximum number of time steps was reached (see MAX_ITER in misc.hpp)
   */
  
  time(&start);
  //time's iteration beginning 
  
  while (T > cur_time  && n < MAX_ITER+1){
    
    if (1 == verif){
      dt1=dt_max;
      
      // save the data in huz_evolution.dat
      if (cur_time >= T_output){
        out->write(h,u,v,z,cur_time);
        T_output+=dt_output;
      }// end if
      
      /* save specific points */
      out_specific_points->save(h,u,v,cur_time);

      //boundary conditions
      boundary(h,u,v,cur_time,NXCELL,NYCELL);
      
      for (int i=1 ; i<NXCELL+1 ; i++){
        for (int j=1 ; j<NYCELL+1 ; j++){
          if (h[i][j]<=HE_CA) {
            h[i][j]=0.;
            u[i][j] = 0.;
            v[i][j] = 0.;
          }
          if (fabs(u[i][j])<=VE_CA){
            u[i][j]=0.;
            q1[i][j]=0.;
          }
          if (fabs(v[i][j])<=VE_CA){
            v[i][j]=0.;
            q2[i][j]=0.;
          }
        } //end for j
      } //end for i
      
      // Reconstruction for order 2
      rec->calcul(h,u,v,z,delzc1,delzc2,delz1,delz2,h1r,u1r,v1r,h1l,u1l,v1l,h2r,u2r,v2r,h2l,u2l,v2l);
      
    }
    else{//if verif==0, ie dt2<dt1
    /*We go back to the values of the previous time */
    for (int i=1 ; i<NXCELL+1 ; i++)
      for (int j=1 ; j<NYCELL+1 ; j++)
        Vin1[i][j]=Vin2[i][j];
    }//end for if verif==1
    
    maincalcflux(CFL_FIX,T, cur_time, dt_max, dt1, dt1);
    
    dt1=min(T-cur_time,dt1);
    
    tx=dt1/DX;
    ty=dt1/DY;
    
    maincalcscheme(h,u,v,q1,q2,hs,us,vs,qs1,qs2,Vin1,cur_time,dt1,n);
    //    dt2=dt_max;
    dt2=dt1;
    
    //boundary conditions
    boundary(hs,us,vs,cur_time+dt1,NXCELL,NYCELL);
    
    for (int i=1 ; i<NXCELL+1 ; i++){
      for (int j=1 ; j<NYCELL+1 ; j++){
        if (hs[i][j]<=HE_CA){
          hs[i][j]=0.;
          us[i][j] = 0.;
          vs[i][j] = 0.;
          qs1[i][j]=0.;
          qs2[i][j]=0.;
        }
        if (fabs(us[i][j])<=VE_CA){
          us[i][j]=0.;
          qs1[i][j]=0.;
        }
        if (fabs(vs[i][j])<=VE_CA){
          vs[i][j]=0.;
          qs2[i][j]=0.;
        }
      } //end for j
    } //end for i
    
    //Reconstruction for order 2 
    rec->calcul(hs,us,vs,z,delzc1,delzc2,delz1,delz2,h1r,u1r,v1r,h1l,u1l,v1l,h2r,u2r,v2r,h2l,u2l,v2l);
    
    //common part
    
    maincalcflux(CFL_FIX,T,cur_time+dt1,dt_max,dt2,dt2);
    
    if (dt2<dt1){
      dt1=dt2;
      tx=dt1/DX;
      ty=dt1/DY;
      verif=0;
    }
    else{
      
      //Added to do calculus at the beginning 
      verif=1;
      maincalcscheme(hs,us,vs,qs1,qs2,hsa,usa,vsa,qsa1,qsa2,Vin1,cur_time,dt1,n+1);
      
      /*the values of height_Vinf_tot and height_of_tot are put to zero
       to compute infiltrated and overland flow volume*/
      height_Vinf_tot=ZERO;
      height_of_tot=ZERO;
      //Heun method (order 2 in time)
      for (int i=1 ; i<NXCELL+1 ; i++){
        for (int j=1 ; j<NYCELL+1 ; j++){
          if (hsa[i][j]<=HE_CA){
            hsa[i][j]=0.;
          }
          if (abs(usa[i][j])<=VE_CA){
            usa[i][j]=0.;
          }
          if (abs(vsa[i][j])<=VE_CA){
            vsa[i][j]=0.;
          }
          tmp = 0.5*(h[i][j]+hsa[i][j]);
          if (tmp>=HE_CA){
            q1[i][j] = 0.5*(h[i][j]*u[i][j]+hsa[i][j]*usa[i][j]);
            u[i][j] = q1[i][j]/tmp;
            q2[i][j] = 0.5*(h[i][j]*v[i][j]+hsa[i][j]*vsa[i][j]);
            v[i][j] = q2[i][j]/tmp;
            h[i][j] = tmp;
          }
          else{
            u[i][j] = 0.;
            q1[i][j] = 0.;
            v[i][j] = 0.;
            q2[i][j] = 0.;
            h[i][j] = 0.;
          }
          Vin_tot[i][j] = (Vin1[i][j] + Vin2[i][j])*0.5;
          
          Vin1[i][j]=Vin_tot[i][j];
          Vin2[i][j]=Vin_tot[i][j];
          
          height_Vinf_tot+= Vin_tot[i][j];
          height_of_tot+=h[i][j];
        } //end for j
      } //end for i
      
      /*Computation of the cumulated infiltrated volume*/
      Vol_inf_tot_cumul=height_Vinf_tot*DX*DY;
      
      /*Computation of the overland flow volume*/
      Vol_of_tot=height_of_tot*DX*DY;
      
      cur_time=cur_time+dt1;
      n=n+1;
      
#ifdef DEBUG
      out->check_vol(cur_time,dt1,Volrain_Tot,Vol_inf_tot_cumul,Vol_of_tot,Total_volume_outflow);
#endif
      
      //Displays the percentage of elapsed time
      cout  << '\r' << '\t' << "[" << int((cur_time/T)*100) << "%] done" ;
      
    }//end for else dt2<dt1
    
  } //end for while : loop in time
  
  //Verifies the reason why the run finished
  if (n>MAX_ITER){
    cerr << "\n order2: WARNING: the computation finished because the maximum number of time steps was reached (see MAX_ITER in misc.hpp)\n";
  }
  
  //Saves the last results
  out->write(h,u,v,z,cur_time);
  
  /* save the last for specific points */
  out_specific_points->save(h,u,v,cur_time);
  
  //to give the computing time
  time(&end);
  timecomputation=difftime(end,start);
  cpu_time = clock();
  
  //to inform about the froude number
  Fr=froude_number(hs,us,vs);
  
  
  // The quantity of water outflow is the sum of flux at the boundary multiply by one cell area, so
  //Outflow volum = (fluxNxcell_cum_T+fluxNycell_cum_T)*DX*DY*n1+(fluxx0_cum_T+fluxy0_cum_T)*DX*DY*n2
  //In this case n1=1 and n2=-1 because we consider the direction of the flow
  // is from left to the right (x=0 to x=Nxcell) and from bottom to the top (y=0 to y=NYCELL)  
  out->result(timecomputation,cpu_time,Volrain_Tot, Vol_inf_tot_cumul, Vol_of_tot,Fr,n,Total_volume_outflow);
  //storage of h u v value in the final time
  out->final(z, h, u,v);
  
}