Traitement_particulier_NS_EC.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_EC.h>
#include <Navier_Stokes_std.h>
#include <Domaine_VF.h>
#include <Terme_Source_Acceleration.h>
#include <Milieu_base.h>
#include <TRUST_Ref.h>
#include <Probleme_base.h>
#include <Schema_Temps_base.h>
#include <sys/stat.h>
 
Implemente_base_sans_constructeur_ni_destructeur(Traitement_particulier_NS_EC,"Traitement_particulier_NS_EC",Traitement_particulier_NS_base);
// XD ec traitement_particulier_base ec INHERITS_BRACE Keyword to print total kinetic energy into the referential linked
// XD_CONT to the domain (keyword Ec). In the case where the domain is moving into a Galilean referential, the keyword
// XD_CONT Ec_dans_repere_fixe will print total kinetic energy in the Galilean referential whereas Ec will print the
// XD_CONT value calculated into the moving referential linked to the domain
// XD attr Ec rien Ec OPT not_set
// XD attr Ec_dans_repere_fixe rien Ec_dans_repere_fixe OPT not_set
// XD attr periode floattant periode OPT periode is the keyword to set the period of printing into the file
// XD_CONT datafile_Ec.son or datafile_Ec_dans_repere_fixe.son.
 
/*! @brief Prints the object to an output stream.
 *
 * @param is an output stream
 * @return the modified output stream
 */
Sortie& Traitement_particulier_NS_EC::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_EC::readOn(Entree& is)
{
  return is;
}
 
Entree& Traitement_particulier_NS_EC::lire(Entree& is)
{
  Motcle accouverte = "{" , accfermee = "}" ;
  //Motcle valec="Ec";
  Motcle motbidon, motlu;
  periode = mon_equation->probleme().schema_temps().pas_temps_min();
  is >> motbidon ;
  if (motbidon == accouverte)
    {
      Motcles les_mots(2);
      les_mots[0] = "Ec";
      les_mots[1] = "Ec_dans_repere_fixe";
      {
        is >> motlu;
        while(motlu != accfermee)
          {
            int rang=les_mots.search(motlu);
            switch(rang)
              {
              case 0 :
              case 1 :
                {
                  if (rang==1)
                    repere_mobile_=1;
                  is >> motlu;
                  if (motlu!="periode")
                    {
                      Cerr<<"We expected the keyword periode after " << les_mots[rang] << "."<<finl;
                      exit();
                    }
                  is >> periode;
                  break;
                }
              default :
                {
                  Cerr << "Error while reading Traitement_particulier_NS_EC";
                  Cerr << "Possible keywords are "<< les_mots <<" { and }" << finl;
                  Cerr << "You read:" << motlu << finl;
                  exit();
                  break;
                }
              }
            is >> motlu;
          }
        is >> motlu;
        if (motlu != accfermee)
          {
            Cerr << "Error while reading Traitement_particulier_NS_EC";
            Cerr << "We expected a }" << finl;
            exit();
          }
      }
    }
  else
    {
      Cerr << "Error while reading Traitement_particulier_NS_EC";
      Cerr << "We expected a {" << finl;
      exit();
    }
  return is;
}
 
void Traitement_particulier_NS_EC::ouvrir_fichier(SFichier& s) const
{
  // Open the file Nom_du_cas_EC.son
  Nom nom_fich(nom_du_cas());
  nom_fich += "_EC";
  if (repere_mobile_)
    nom_fich+="_dans_repere_fixe";
  nom_fich+=".son";
 
  const Probleme_base& pb=mon_equation->probleme();
  struct stat f;
  if (stat(nom_fich,&f) || (nb_bip==0 && !pb.reprise_effectuee()))
    s.ouvrir(nom_fich);
  else
    s.ouvrir(nom_fich,ios::app);
  s.setf(ios::scientific);
  s.precision(8);
}
 
void Traitement_particulier_NS_EC::preparer_calcul_particulier()
{
  tinit = mon_equation->inconnue().temps();
  nb_bip=0;
  double Ec;
  calculer_Ec(Ec);
  if(je_suis_maitre())
    {
      SFichier le_fichier;
      ouvrir_fichier(le_fichier);
      // Write first line
      le_fichier<<"# Temps        Energie_cinetique_totale"<<finl;
      le_fichier<<tinit<<" "<<Ec<<finl;
    }
}
void Traitement_particulier_NS_EC::post_traitement_particulier()
{
  double temps = mon_equation->inconnue().temps();
  double nb = floor((temps - tinit) / periode);
  if (nb > nb_bip + 0.5)
    {
      nb_bip=nb;
      double Ec;
      calculer_Ec(Ec);
      if(je_suis_maitre())
        {
          SFichier le_fichier;
          ouvrir_fichier(le_fichier);
          le_fichier<<temps<<" "<<Ec<<finl;
        }
    }
}
 
/*! @brief Helper function used in calculer_Ec.
 *
 * Computes the sum of 0.5*v^2*rho*volumes_entrelaces.
 *
 */
static double trait_part_calculer_ec_faces(const int         face_debut,
                                           const int         nb_faces,
                                           const int                frontiere,
                                           const DoubleTab&         vitesse,
                                           const DoubleVect&         volumes_entrelaces,
                                           const DoubleTab&         xv,
                                           const DoubleTab&         masse_volumique,
                                           const ArrOfDouble&         translation,
                                           const ArrOfDouble&         rotation,
                                           const int         repere_mobile_,
                                           const ArrOfInt&        faces_doubles
                                          )
{
  const int face_fin = face_debut + nb_faces;
  double ec = 0.;
  double rho = 0.;
  const int nb_dim_1 = (vitesse.line_size() == 1);
  const int dim      = Objet_U::dimension;
  ArrOfDouble ve(Objet_U::dimension);
  for (int face = face_debut; face < face_fin; face++)
    {
      // Computation of the entrainment velocity
      if (repere_mobile_)
        {
          ve[0]=translation[0];
          ve[1]=translation[1];
          if (Objet_U::dimension==3)
            {
              ve[2]=translation[2];
              ve[0]+=rotation[1]*xv(face,2)-rotation[2]*xv(face,1);
              ve[1]+=rotation[2]*xv(face,0)-rotation[0]*xv(face,2);
              ve[2]+=rotation[0]*xv(face,1)-rotation[1]*xv(face,0);
            }
        }
      else
        ve=0;
 
      double v2;
      double volume;
      if (nb_dim_1)
        {
          // One velocity component at the face (VDF)
          const double v = vitesse(face);
          if (repere_mobile_)
            {
              Cerr << "The computation of kinetic energy in a fixed frame" <<finl;
              Cerr << "is not implemented in VDF." << finl;
              Process::exit(); // Indeed, a design issue: the orientation of VDF faces would be needed
            }
          v2 = v * v;
          // In VDF, at boundaries, only half the staggered volume is taken
          volume = (frontiere ? 0.5 : 1) * volumes_entrelaces(face);
        }
      else
        {
          // Two or three components (VEFP1B)
          v2 = 0.;
          for (int i = 0; i < dim; i++)
            {
              const double v_i = vitesse(face, i);
              v2 += (v_i + ve[i]) * (v_i + ve[i]);
            }
          // In VEF, this is incorrect; extended volumes should be used:
          volume = volumes_entrelaces(face);
        }
      const int k = (masse_volumique.dimension(0)==1) ? 0 : face;
      rho = masse_volumique(k, 0);
      double contribution = (faces_doubles[face]==1) ? 0.5 : 1 ;
      ec += contribution * 0.5 * v2 * volume * rho;
    }
  return ec;
}
 
/*! @brief Same computation method for VDF and VEF.
 *
 *   If no mass-density field has been associated, computes
 *    INTEGRAL of 1/2 * u^2 over the domain.
 *   If a mass-density field has been associated, it must be
 *   of type "champ aux faces". Computes INTEGRAL of 1/2 * rho * u^2.
 *   rho is a P0 field at elements, u is P0 on the staggered volumes.
 *  Return value:
 *   The sum over all processors of the integral of 1/2*rho*u*u.
 *
 */
void Traitement_particulier_NS_EC::calculer_Ec(double& energie_cinetique)
{
 
  const Domaine_dis_base& zdisbase        = mon_equation->inconnue().domaine_dis_base();
  const Domaine_VF&    domaine_VF            = ref_cast(Domaine_VF, zdisbase);
  const DoubleVect& volumes_entrelaces = domaine_VF.volumes_entrelaces();
  const DoubleTab&  xv                 = domaine_VF.xv();
  const DoubleTab&  vitesse            = mon_equation->inconnue().valeurs();
  OBS_PTR(ArrOfDouble) translation(xv);
  OBS_PTR(ArrOfDouble) rotation(xv);
  OBS_PTR(DoubleTab) rho(xv);
  DoubleVect rotation_nulle(dimension);
  rotation_nulle=0;
  if (repere_mobile_)
    {
      int ok=0;
      // Check for the existence of an Acceleration source term in Navier Stokes
      const Sources& les_sources=mon_equation->sources();
 
      for (const auto& itr : les_sources)
        {
          if (sub_type(Terme_Source_Acceleration,itr.valeur()))
            {
              const Terme_Source_Acceleration& terme_source_acceleration=ref_cast(Terme_Source_Acceleration,itr.valeur());
              // Check that the translation velocities of the moving
              // frame are properly defined and associate the translation and
              // optionally rotation arrays
              if (terme_source_acceleration.has_champ_vitesse())
                {
                  translation=terme_source_acceleration.champ_vitesse().valeurs();
                  if (terme_source_acceleration.has_omega())
                    rotation=terme_source_acceleration.omega().valeurs();
                  else
                    rotation=rotation_nulle;
                  ok=1;
                }
            }
        }
      if (!ok)
        {
          Cerr << "You cannot compute kinetic energy in a fixed frame" << finl;
          Cerr << "unless the calculation frame is moving, i.e. you have" << finl;
          Cerr << "defined an acceleration source term in the Navier Stokes equation." << finl;
          Cerr << "Or the definition of the velocity of the" << finl;
          Cerr << "moving frame in the fixed frame is missing in this source term. Modify your data file." << finl;
          exit();
        }
    }
  double ec = 0.;
 
  if (has_champ_masse_volumique())
    {
      const Champ_base& champ_rho = get_champ_masse_volumique();
      rho = champ_rho.valeurs();
      if (rho->dimension(0) != domaine_VF.nb_faces() || rho->line_size() != 1)
        {
          Cerr << "Error in Traitement_particulier_NS_EC::calculer_Ec" << finl;
          Cerr << "the density field is not a scalar field at faces" << finl;
          Process::exit();
        }
    }
  else
    {
      rho = mon_equation->milieu().masse_volumique().valeurs();
    }
  const int nb_front = domaine_VF.nb_front_Cl();
  const ArrOfInt& faces_doubles = domaine_VF.faces_doubles();
  // Boundary faces
  for (int i = 0; i < nb_front; i++)
    {
      const Frontiere& fr = domaine_VF.front_VF(i).frontiere();
      const int debut    = fr.num_premiere_face();
      const int nb_faces = fr.nb_faces();
      ec += trait_part_calculer_ec_faces(debut, nb_faces, 1,
                                         vitesse, volumes_entrelaces, xv, rho, translation, rotation, repere_mobile_, faces_doubles);
    }
  // Interior faces of the domain (including joint faces)
  {
    const int debut    = domaine_VF.premiere_face_int();
    const int nb_faces = domaine_VF.nb_faces_internes();
    ec += trait_part_calculer_ec_faces(debut, nb_faces, 0,
                                       vitesse, volumes_entrelaces, xv, rho, translation, rotation, repere_mobile_, faces_doubles);
  }
  ec = Process::mp_sum(ec);
  energie_cinetique = ec;
}