Traitement_particulier_NS_Profils.cpp

/****************************************************************************
* Copyright (c) 2024, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Traitement_particulier_NS_Profils.h>
#include <LecFicDistribueBin.h>
#include <EcrFicCollecteBin.h>
#include <Modele_turbulence_hyd_base.h>
#include <Navier_Stokes_std.h>
 
Implemente_base(Traitement_particulier_NS_Profils,"Traitement_particulier_NS_Profils",Traitement_particulier_NS_base);
 
 
/*! @brief Prints the object to an output stream.
 *
 * @param is an output stream
 * @return the modified output stream
 */
Sortie& Traitement_particulier_NS_Profils::printOn(Sortie& is) const
{
  return is;
}
 
 
/*! @brief Reads the object from an input stream.
 *
 * @param is an input stream
 * @return the modified input stream
 */
Entree& Traitement_particulier_NS_Profils::readOn(Entree& is)
{
  return is;
}
 
Entree& Traitement_particulier_NS_Profils::lire(Entree& is)
{
  // Initialisation
  IntVect Verif(8);
  Verif=0;
 
  // END Initialisation
  Motcle accouverte = "{" , accfermee = "}" ;
  Motcle motbidon, motlu;
  is >> motbidon ;
  if (motbidon == accouverte)
    {
      Motcles les_mots(8);
      les_mots[0] = "const_profile";
      les_mots[1] = "num_probes";
      les_mots[2] = "dt_output";
      les_mots[3] = "stats";
      les_mots[4] = "profil_nu_t";
      les_mots[5] = "reprise";
      les_mots[6] = "reprise_nu_t";
      les_mots[7] = "homog_dir";
      is >> motlu;
      while(motlu != accfermee)
        {
          int rang=les_mots.search(motlu);
          switch(rang)
            {
            case 0 :
              {
                // For const_profile Constant Profile value : X, Y, Z
                is >> motlu;      // Direction of constant profile
                if((motlu=="x")||(motlu=="X")) dir_profil = 0;
                else if((motlu=="y")||(motlu=="Y")) dir_profil = 1;
                else if((motlu=="z")||(motlu=="Z")) dir_profil = 2;
                else
                  {
                    Cerr << "Wrong value for direction of profile !" << finl;
                    exit();
                  }
                Cerr << "Profiles at constant " << motlu << finl;
                break;
              }
            case 1 :
              {
                // Number of probes and reading of the values of the postions
                is >> n_probes;
                if(n_probes<1)
                  {
                    Cerr << "Error in Traitement_particulier_NS_Profils::lire" << finl;
                    Cerr << "Wrong value for number of probes (" << n_probes << ") :" << finl;
                    Cerr << "There must be at least one profile !" << finl;
                    exit();
                  }
                positions.resize(n_probes);
                for(int i=0; i<n_probes; i++)
                  {
                    is >> positions(i);
                    Cerr << "  Probe_" << i+1 << " = " << positions(i) << finl;
                  }
                break;
              }
            case 2 :
              {
                // For u_inst
                is >> dt_post_inst;      // time interval for printing spatial averages
                // initialisation for spatial stats computation
                oui_u_inst = 1;
                Cerr << "Spatial average parameters read..." << finl;
                Cerr << "dt_post_inst : " << dt_post_inst << finl;
                break;
              }
            case 3 :
              {
                // For Time average
                is >> temps_deb;     // Initial time for time average statistics.
                is >> temps_fin;     // Final time for time average statistics.
                is >> dt_post_stat;  // Output intervalfor time average statistics.
                oui_stat=1;          // =1 : computing temporal stats
                Cerr << "Reading time statitics parameters..." << finl;
                Cerr << "Initial time : " << temps_deb << " End time : " << temps_fin << finl;
                Cerr << "Time interval for output : " << dt_post_stat << finl;
                // Check whether values specific to spatial stats computations have been provided
                break;
              }
            case 4 :
              {
                // For eddy viscosity profile
                oui_profil_nu_t = 1;
                Cerr << "User asked for statistics on eddy-viscosity..." << finl;
                const RefObjU& modele_turbulence = mon_equation->get_modele(TURBULENCE);
                const Modele_turbulence_hyd_base& mod_turb = ref_cast(Modele_turbulence_hyd_base,modele_turbulence.valeur());
                if(sub_type(Modele_turbulence_hyd_base,mod_turb))
                  Cerr << "Statistics for eddy-viscosity : OK !" << finl;
                else
                  {
                    Cerr << "User asked for eddy-viscosity profile," << finl;
                    Cerr << "but is not defined in laminar case !" << finl;
                    Cerr << "Please check the problem type..." << finl;
                    exit();
                  }
                break;
              }
            case 5 :
              {
                oui_repr = 1;
                is  >> fich_repr_u ;  // name of the stats restart file
                Cerr << "The time statistics file is : " << fich_repr_u << finl;
                break;
              }
            case 6 :
              {
                // For eddy viscosity profile Reprise
                oui_repr_nu_t = 1;
                is  >> fich_repr_nu_t ;  // name of the nut stats restart file
                Cerr << "Continuing statistics on eddy-viscosity : OK!" << finl;
                break;
              }
            case 7 :
              {
                // For homogeneous direction specification : X, Y, Z
                is >> motlu;      // Direction of homogeneity
                if((motlu=="x")||(motlu=="X")) homo_dir = 0;
                else if((motlu=="y")||(motlu=="Y")) homo_dir = 1;
                else if((motlu=="z")||(motlu=="Z")) homo_dir = 2;
                else
                  {
                    Cerr << "Wrong value for homogeneous direction !" << finl;
                    exit();
                  }
                Cerr << "Homogeneous direction : " << motlu << finl;
                break;
              }
            default :
              {
                lire(motlu,is);
                break;
              }
            }
          is >> motlu;
        }
      is >> motlu;
      if (motlu != accfermee)
        {
          Cerr << "Error while reading in Traitement_particulier_NS_Profils::lire" << finl;
          Cerr << "We were expecting a }" << finl;
          exit();
        }
    }
  else
    {
      Cerr << "Error while reading in Traitement_particulier_NS_Profils" << finl;
      Cerr << "We were expecting a {" << finl;
      exit();
    }
 
  if ((oui_u_inst!=1)&&(oui_stat!=0))
    {
      Cerr << "User asked for time statistics :" << finl;
      Cerr << "To do so, you also need to specify the parameters for space statistics," << finl;
      Cerr << "since it is used for the time stats.." << finl;
      exit();
    }
  if ((oui_profil_nu_t!=1)&&(oui_repr_nu_t!=0))
    {
      Cerr << "User asked for time statistics on nu_t :" << finl;
      Cerr << "to do so, you also need to request the computation of space statistics," << finl;
      Cerr << "since it is necessary for computing the temporal stats." << finl;
      exit();
    }
  return is;
}
 
// Entree& Traitement_particulier_NS_canal::lire(const Motcle& motlu, Entree& is)
// {
//   return is;
// }
 
void Traitement_particulier_NS_Profils::preparer_calcul_particulier()
{
  if ((oui_u_inst != 0)||(oui_profil_nu_t != 0))
    init_calcul_moyenne();
  if (oui_stat != 0)
    init_calcul_stats();
}
 
void Traitement_particulier_NS_Profils::init_calcul_stats()
{
  if (oui_repr!=1) // this is not a stats restart, so no reading from a file
    {
      u_moy_temp_x.resize(n_probes,Nap);
      u_moy_temp_x=0;
      u_moy_temp_y.resize(n_probes,Nap);
      u_moy_temp_y=0;
      u_moy_temp_z.resize(n_probes,Nap);
      u_moy_temp_z=0;
 
      uv_moy_temp.resize(n_probes,Nap);
      uv_moy_temp=0; //  uv_moy time average
 
      u2_moy_temp.resize(n_probes,Nap);
      u2_moy_temp=0;
      v2_moy_temp.resize(n_probes,Nap);
      v2_moy_temp=0;
      w2_moy_temp.resize(n_probes,Nap);
      w2_moy_temp=0;
 
      u3_moy_temp.resize(n_probes,Nap);
      u3_moy_temp=0;
      v3_moy_temp.resize(n_probes,Nap);
      v3_moy_temp=0;
      w3_moy_temp.resize(n_probes,Nap);
      w3_moy_temp=0;
 
      u4_moy_temp.resize(n_probes,Nap);
      u4_moy_temp=0;
      v4_moy_temp.resize(n_probes,Nap);
      v4_moy_temp=0;
      w4_moy_temp.resize(n_probes,Nap);
      w4_moy_temp=0;
 
      nu_t_temp.resize(n_probes,Nap);
      nu_t_temp=0;
    }
}
 
void Traitement_particulier_NS_Profils::reprendre_stat()
{
  double tps = mon_equation->inconnue().temps();
  int j,num_prof;
  double ti,ti2;
  Nom temps;
  if (oui_stat == 1)
    {
      ifstream ficu0(fich_repr_u);
      if (ficu0)
        {
          LecFicDistribueBin ficu (fich_repr_u);
          ficu >> temps;
          u_moy_temp_x.resize(n_probes,Nap);
          u_moy_temp_y.resize(n_probes,Nap);
          u_moy_temp_z.resize(n_probes,Nap);
 
          uv_moy_temp.resize(n_probes,Nap);
          u2_moy_temp.resize(n_probes,Nap);
          v2_moy_temp.resize(n_probes,Nap);
          w2_moy_temp.resize(n_probes,Nap);
 
          u3_moy_temp.resize(n_probes,Nap);
          v3_moy_temp.resize(n_probes,Nap);
          w3_moy_temp.resize(n_probes,Nap);
 
          u4_moy_temp.resize(n_probes,Nap);
          v4_moy_temp.resize(n_probes,Nap);
          w4_moy_temp.resize(n_probes,Nap);
 
          for(num_prof=0; num_prof<n_probes; num_prof++)
            {
              ficu >> Nxy(num_prof) >> Nyy(num_prof) >> Nzy(num_prof) >> Nuv(num_prof);
              for (j=0; j<Nxy(num_prof); j++)
                {
                  ficu >> u_moy_temp_x(num_prof,j) ;
                  ficu >> u2_moy_temp(num_prof,j);
                }
              for (j=0; j<Nyy(num_prof); j++)
                {
                  ficu >> u_moy_temp_y(num_prof,j) ;
                  ficu >> v2_moy_temp(num_prof,j);
                }
              for (j=0; j<Nzy(num_prof); j++)
                {
                  ficu >> u_moy_temp_z(num_prof,j);
                  ficu >> w2_moy_temp(num_prof,j);
                }
              for (j=0; j<Nuv(num_prof); j++)
                ficu >> uv_moy_temp(num_prof,j) ;
            }
          //END reading "n_probes" profiles
          ficu >> ti ;
          Cerr << "ti=" << ti << finl;
          Nom chc_ti = Nom(ti), chc_temps_deb = Nom(temps_deb);
          Cerr << "chc_ti=" << chc_ti << "   chc_temps_deb=" << chc_temps_deb << finl;
          if (chc_ti!=chc_temps_deb)
            {
              if (ti > temps_deb)
                {
                  Cerr << "problem resuming velocity stats!!" << finl;
                  Cerr << "The requested start time for stats " << temps_deb  << finl;
                  Cerr << "is less than the start time of the saved stats!" << ti << finl;
                  exit();
                }
              else
                {
                  if (temps_deb>=tps)
                    {
                      Cerr << "Restarting velocity stats computation from t=" << temps_deb << "s" << finl;
                      u_moy_temp_x = 0.;
                      u_moy_temp_y = 0.;
                      u_moy_temp_z = 0.;
                      uv_moy_temp = 0.;
                      u2_moy_temp = 0.;
                      v2_moy_temp = 0.;
                      w2_moy_temp = 0.;
                      oui_repr = 0;
                    }
                  else
                    {
                      Cerr << "The time at which you want to start stats has already been passed!!" << finl;
                      exit();
                    }
                }
            }
          else
            {
              Cerr << "Continuing velocity stats computation, started at t=" <<  ti << "s" << finl;
              temps_deb = ti;
            }
        }
      else
        {
          if (tps<=temps_deb)
            {
              Cerr << "Stats computation has not started yet" << finl;
              oui_repr = 0;
            }
          else
            {
              Cerr << "You must provide the file to resume the stats!!" << finl;
              exit();
            }
        }
 
      if ((oui_profil_nu_t!=0)&&(oui_repr_nu_t!=0))
        {
          ifstream ficu1(fich_repr_nu_t);
          if (ficu1)
            {
              LecFicDistribueBin ficu (fich_repr_nu_t);
              ficu >> temps;
              nu_t_temp.resize(n_probes,Nap);
              for(num_prof=0; num_prof<n_probes; num_prof++)
                {
                  ficu >>  Nuv(num_prof);
                  for (j=0; j<Nuv(num_prof); j++)
                    ficu >> nu_t_temp(num_prof,j) ;
                }
              ficu >> ti2 ;
              if (ti2!=ti)
                {
                  Cerr << "Problem : the nut stats and velocity stats are misaligned!!" << finl;
                  exit();
                }
              Cerr << "ti2=" << ti2 << finl;
              Nom chc_ti2 = Nom(ti2), chc_temps_deb = Nom(temps_deb);
              Cerr << "chc_ti2=" << chc_ti2 << "   chc_temps_deb=" << chc_temps_deb << finl;
              if (chc_ti2!=chc_temps_deb)
                {
                  if (ti2 > temps_deb)
                    {
                      Cerr << "problem resuming nu_t stats!!" << finl;
                      Cerr << "The requested start time for stats " << temps_deb  << finl;
                      Cerr << "is less than the start time of the saved stats!" << ti << finl;
                      exit();
                    }
                  else
                    {
                      if (temps_deb>=tps)
                        {
                          Cerr << "Restarting nu_t stats computation from t=" << temps_deb << "s" << finl;
                          nu_t_temp = 0.;
                          oui_repr_nu_t = 0;
                        }
                      else
                        {
                          Cerr << "The time at which you want to start stats has already been passed!" << finl;
                          exit();
                        }
                    }
                }
              else
                {
                  Cerr << "Continuing nu_t stats computation, started at t=" <<  ti << "s" << finl;
                  temps_deb = ti2;
                }
            }
          else
            {
              if (tps<=temps_deb)
                {
                  Cerr << "Stats computation has not started yet" << finl;
                  oui_repr_nu_t = 0;
                }
              else
                {
                  Cerr << "You must provide the file to resume the stats!" << finl;
                  exit();
                }
            }
        }
    }
  else
    {
      Cerr << "No stats computation started!" << finl;
      Cerr << "No resumption of u_moy_temp, u_et_temp and teta_moy_temp" << finl;
    }
}
 
void Traitement_particulier_NS_Profils::sauver_stat() const
{
  double tps = mon_equation->inconnue().temps();
  if (  (oui_stat == 1)&&(tps>=temps_deb)&&(tps<=temps_fin) )
    {
      Cerr << "In Traitement_particulier_NS_canal::sauver_stat!!" << finl;
      // Save the sum (without dividing by dt)
      int j;
      Nom temps = Nom(tps);
      Nom fich_sauv_u ="u_moy_temp_";
      fich_sauv_u+=temps;
      fich_sauv_u+=".sauv";
 
      // Save u_moy!!
      EcrFicCollecteBin ficu (fich_sauv_u);
      //      EcrFicCollecte ficu (fich_sauv_u);
      //      EcrFicCollecteBin ficu ("u_moy_temp.sauv");
      ficu << temps << finl;
      int num_prof;
      for(num_prof=0; num_prof<n_probes; num_prof++)
        {
          ficu << Nxy(num_prof) << Nyy(num_prof) << Nzy(num_prof) << Nuv(num_prof) << finl;
          for (j=0; j<Nxy(num_prof); j++)
            {
              ficu << u_moy_temp_x(num_prof,j) << finl;
              ficu << u2_moy_temp(num_prof,j) << finl;
            }
          for (j=0; j<Nyy(num_prof); j++)
            {
              ficu << u_moy_temp_y(num_prof,j) << finl;
              ficu << v2_moy_temp(num_prof,j) << finl;
            }
          for (j=0; j<Nzy(num_prof); j++)
            {
              ficu << u_moy_temp_z(num_prof,j) << finl;
              ficu << w2_moy_temp(num_prof,j) << finl;
            }
          for (j=0; j<Nuv(num_prof); j++)
            ficu << uv_moy_temp(num_prof,j) << finl;
        }
      ficu << temps_deb << finl;
      Cerr << "Saving u_moy_temp at t=" << tps << finl;
      Cerr << "Start of stats computation at t=" << temps_deb << " to note for stats resumption!" << finl;
 
      if (oui_profil_nu_t!=0)
        {
          Nom fich_sauv_nut ="nut_temp_";
          fich_sauv_nut+=temps;
          fich_sauv_nut+=".sauv";
          // Save nut!!
          EcrFicCollecteBin ficu2 (fich_sauv_nut);
          ficu2 << temps << finl;
          for(num_prof=0; num_prof<n_probes; num_prof++)
            {
              ficu2 << Nuv(num_prof) << finl;
              for (j=0; j<Nuv(num_prof); j++)
                ficu2 << nu_t_temp(num_prof,j) << finl;
            }
          ficu2 << temps_deb << finl;
          Cerr << "Saving nu_t_temp at t=" << tps << finl;
          Cerr << "Start of stats computation at t=" << temps_deb << " to note for stats resumption!" << finl;
        }
    }
}