Terme_Source_Acceleration.cpp

/****************************************************************************
* Copyright (c) 2026, 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 <Terme_Source_Acceleration.h>
 
#include <Probleme_base.h>
#include <Navier_Stokes_std.h>
#include <Discretisation_base.h>
#include <Domaine_VF.h>
#include <Schema_Temps_base.h>
 
Implemente_base_sans_constructeur(Terme_Source_Acceleration,"Terme_Source_Acceleration",Source_base);
// XD acceleration source_base acceleration BRACE Momentum source term to take in account the forces due to rotation or
// XD_CONT translation of a non Galilean referential R\' (centre 0\') into the Galilean referential R (centre 0).
// XD attr vitesse field_base vitesse OPT Keyword for the velocity of the referential R\' into the R referential
// XD_CONT (dOO\'/dt term [m.s-1]). The velocity is mandatory when you want to print the total cinetic energy into the
// XD_CONT non-mobile Galilean referential R (see Ec_dans_repere_fixe keyword).
// XD attr acceleration field_base acceleration OPT Keyword for the acceleration of the referential R\' into the R
// XD_CONT referential (d2OO\'/dt2 term [m.s-2]). field_base is a time dependant field (eg: Champ_Fonc_t).
// XD attr omega field_base omega OPT Keyword for a rotation of the referential R\' into the R referential [rad.s-1].
// XD_CONT field_base is a 3D time dependant field specified for example by a Champ_Fonc_t keyword. The time_field field
// XD_CONT should have 3 components even in 2D (In 2D: 0 0 omega).
// XD attr domegadt field_base domegadt OPT Keyword to define the time derivative of the previous rotation [rad.s-2].
// XD_CONT Should be zero if the rotation is constant. The time_field field should have 3 components even in 2D (In 2D:
// XD_CONT 0 0 domegadt).
// XD attr centre_rotation field_base centre_rotation OPT Keyword to specify the centre of rotation (expressed in R\'
// XD_CONT coordinates) of R\' into R (if the domain rotates with the R\' referential, the centre of rotation is
// XD_CONT 0\'=(0,0,0)). The time_field should have 2 or 3 components according the dimension 2 or 3.
// XD attr option chaine(into=["terme_complet","coriolis_seul","entrainement_seul"]) option OPT Keyword to specify the
// XD_CONT kind of calculation: terme_complet (default option) will calculate both the Coriolis and centrifugal forces,
// XD_CONT coriolis_seul will calculate the first one only, entrainement_seul will calculate the second one only.
 
 
Terme_Source_Acceleration::Terme_Source_Acceleration()
{
  /*
    Noms& nom=champs_compris_.liste_noms_compris();
    nom.dimensionner(1);
    nom[0]="ACCELERATION_TERME_SOURCE";
  */
}
 
Sortie& Terme_Source_Acceleration::printOn(Sortie& s ) const
{
  Cerr << "Terme_Source_Acceleration::printOn : invalid call" << finl;
  assert(0);
  exit();
  return s;
}
 
/*! @brief Call to Terme_Source_Acceleration::lire_data
 *
 */
Entree& Terme_Source_Acceleration::readOn(Entree& s )
{
  Cerr << "Terme_Source_Acceleration::readOn : invalid call" << finl;
  assert(0);
  exit();
  return s;
}
 
/*! @brief Method called by readOn of derived classes Terme_Source_Acceleration_VDF_Face, .
 *
 * ..
 *
 *  Expected format:
 *   {
 *     [ omega           A_Uniform_Field
 *       domegadt        A_Uniform_Field
 *       centre_rotation A_Uniform_Field
 *       option          TERME_COMPLET | CORIOLIS_SEUL | ENTRAINEMENT_SEUL ]
 *     [ acceleration    A_Uniform_Field ]
 *   }
 *  A_Uniform_Field is typically:
 *       Champ_Uniforme 3 x y z
 *     ou
 *       Champ_Fonc_t DOM 3 fx(t) fy(t) fz(t)
 *
 */
 
void Terme_Source_Acceleration::lire_data(Entree& s)
{
  if (je_suis_maitre())
    Cerr << "Terme_Source_Acceleration::lire_data" << finl;
 
  Motcles les_mots(7);
  les_mots[0] = "omega";
  les_mots[1] = "domegadt";
  les_mots[2] = "centre_rotation";
  les_mots[3] = "option";
  les_mots[4] = "acceleration";
  les_mots[5] = "}";
  les_mots[6] = "vitesse";
 
  Motcle motlu;
  s >> motlu;
  if (motlu != "{")
    {
      Cerr << "Error in Terme_Source_Acceleration::lire_data" << finl;
      Cerr << " Expected {" << finl;
      assert(0);
      exit();
    }
  int fini = 0;
  while (!fini)
    {
      s >> motlu;
      const int i_mot = les_mots.search(motlu);
      if (je_suis_maitre() && i_mot < 5)
        Cerr << "Reading " << motlu << " : ";
      switch(i_mot)
        {
        case 0:
          s >> omega_;
          if (je_suis_maitre())
            Cerr << omega_.que_suis_je() << finl;
          break;
        case 1:
          s >> domegadt_;
          if (je_suis_maitre())
            Cerr << domegadt_.que_suis_je() << finl;
          break;
        case 2:
          s >> centre_rotation_;
          if (je_suis_maitre())
            Cerr << centre_rotation_.que_suis_je() << finl;
          break;
        case 3:
          {
            Motcles les_options(3);
            les_options[0] = "terme_complet";
            les_options[1] = "coriolis_seul";
            les_options[2] = "entrainement_seul";
            s >> motlu;
            const int i = les_options.search(motlu);
            switch (i)
              {
              case 0:
                option_ = TERME_COMPLET;
                break;
              case 1:
                option_ = CORIOLIS_SEUL;
                break;
              case 2:
                option_ = ENTRAINEMENT_SEUL;
                break;
              default:
                Cerr << "Error, valid options are: " << finl;
                Cerr << les_options << finl;
                assert(0);
                exit();
              }
            if (je_suis_maitre())
              Cerr << motlu << finl;
            break;
          }
        case 4:
          s >> champ_acceleration_;
          if (je_suis_maitre())
            Cerr << champ_acceleration_.que_suis_je() << finl;
          break;
        case 5:
          fini = 1;
          break;
        case 6:
          s >> champ_vitesse_;
          if (je_suis_maitre())
            Cerr << champ_vitesse_.que_suis_je() << finl;
          break;
        default:
          Cerr << "Error in Terme_Source_Acceleration::lire_data" << finl;
          Cerr << " Keyword not recognized: " << motlu << finl;
          Cerr << " Recognized keywords are: " << les_mots << finl;
          assert(0);
          exit();
        }
    }
 
  // Check consistency:
  int n = 0;
  if (omega_) n++;
  if (domegadt_) n++;
  if (centre_rotation_) n++;
  if (n != 0 && n != 3)
    {
      Cerr << "Error in Terme_Source_Acceleration::lire_data" << finl;
      Cerr << " If OMEGA or DOMEGADT or CENTRE_GRAVITE is given" << finl;
      Cerr << " then all three must be given" << finl;
      assert(0);
      exit();
    }
  if (n > 0)
    {
      if (omega_->valeurs().dimension(0) != 1
          || domegadt_->valeurs().dimension(0) != 1
          || omega_->valeurs().dimension(1) != 3
          || domegadt_->valeurs().dimension(1) != 3)
        {
          Cerr << "Error in Terme_Source_Acceleration::lire_data" << finl;
          Cerr << " the OMEGA and DOMEGADT fields must be uniform fields" << finl;
          Cerr << " with three components (vector aligned on Z in 2D)" << finl;
          assert(0);
          exit();
        }
    }
  if (champ_acceleration_)
    {
      if (champ_acceleration_->valeurs().dimension(0) != 1)
        {
          Cerr << "Error in Terme_Source_Acceleration::lire_data" << finl;
          Cerr << " The ACCELERATION field must be a uniform field" << finl;
          assert(0);
          exit();
        }
    }
 
  // Discretization of la_source and terme_source_post
  // This must be done here so that the fields are available at the time
  // post-processing is read (call to a_pour_Champ_Fonc)
  {
    const Equation_base& eqn = equation();
    if (!sub_type(Navier_Stokes_std, eqn))
      {
        Cerr << "Error in Terme_Source_Acceleration::lire_data\n"
             << " This source term can only be inserted into an equation\n"
             << " derived from Navier_Stokes_std" << finl;
        assert(0);
        exit();
      }
    const Discretisation_base& dis   = eqn.discretisation();
    const Domaine_dis_base&        domaine  = eqn.domaine_dis();
    const double                temps = eqn.schema_temps().temps_courant();
    dis.discretiser_champ("vitesse", domaine, "acceleration_terme_source", "kg/ms^2",
                          Objet_U::dimension, /* composantes */
                          temps,
                          terme_source_post_);
    champs_compris_.ajoute_champ(terme_source_post_);
  }
}
 
void Terme_Source_Acceleration::associer_pb(const Probleme_base& pb)
{
}
 
const Champ_Fonc_base& Terme_Source_Acceleration::get_terme_source_post() const
{
  return terme_source_post_.valeur();
}
 
Champ_Fonc_base& Terme_Source_Acceleration::get_set_terme_source_post() const
{
  // terme_source_post_ is mutable, so it can be returned as "non-const"
  return terme_source_post_.valeur();
}
 
/*! @brief Computes the value of the la_source field at faces based on - calculer_vitesse_faces()
 *
 *   - champ_acceleration_
 *   - omega_
 *   - domegadt_
 *   - centre_rotation_
 *
 * @param acceleration_aux_faces a field discretized at elements of the domaine_VF in which the result of the computation is stored. Virtual space up to date.
 */
const DoubleTab&
Terme_Source_Acceleration::calculer_la_source(DoubleTab& acceleration_aux_faces) const
{
  const Domaine_VF&    domaine_VF = ref_cast(Domaine_VF, equation().domaine_dis());;
 
  const DoubleTab& centre_faces = domaine_VF.xv();
 
  double a[3] = {0., 0., 0.};  // champ acceleration
  double w[3] = {0., 0., 0.};  // omega
  double dw[3] = {0., 0., 0.}; // domegadt
  double c[3] = {0., 0., 0.};  // centre_gravite
  const int dim = Objet_U::dimension;
 
  // ****************************************************************
  // INITIALIZATION OF a, w, dw, c
  int rotation_solide = 0;
  {
    int j;
 
    if (champ_acceleration_)
      {
        const DoubleTab& a_ = champ_acceleration_->valeurs();
        for (j = 0; j < dim; j++)
          a[j] = a_(0, j);
      }
 
    if (omega_)
      {
        // The user has specified a solid rotation motion
        const DoubleTab& champ_w  = omega_          ->valeurs();
        const DoubleTab& champ_dw = domegadt_       ->valeurs();
        const DoubleTab& champ_c  = centre_rotation_->valeurs();
        // Note: in 2D, the rotation vector is typically oriented
        // in the Z direction: so loop up to 3 in all cases:
        for (j = 0; j < 3; j++) w[j]  = champ_w(0,j);
        for (j = 0; j < 3; j++) dw[j] = champ_dw(0,j);
        for (j = 0; j < dim; j++) c[j]  = champ_c(0,j);
        rotation_solide = 1;
      }
  }
 
  // If we have a solid rotation source term, we need the three
  // components of the fluid velocity at the faces:
 
  // A storage space for the velocity if it needs to be
  // computed (in VDF, the three components are computed at each face)
  DoubleTab vitesse_faces_stockage;
  // Computation of the velocity field at the faces (from the unknown
  // of the N.S. equation)
  const DoubleTab& vitesse_faces =
    (rotation_solide)
    ? calculer_vitesse_faces(vitesse_faces_stockage)
    : vitesse_faces_stockage;
 
  // **************************************************
  // LOOP OVER FACES
  int i_face;
  const int nb_faces = acceleration_aux_faces.dimension(0);
 
  // Fluid velocity vector in the current element
  //  (note: Z component initialized for 2D)
  double v[3] = {0., 0., 0.};
  // Vector (centre_face - centre_rotation)
  //  (same)
  double xgr[3] = {0., 0., 0.};
 
  for (i_face = 0; i_face < nb_faces; i_face++)
    {
 
      int j;
      double src[3];
      src[0] = -a[0]; // Terme d'acceleration
      src[1] = -a[1];
      src[2] = -a[2]; // =0. en 2D
 
      // Formules identiques en 2D et en 3D :
      //  en 2D, les composantes z de v et xgr sont nulles.
      if (rotation_solide)
        {
          for (j = 0; j < dim; j++) // en 2D, v[3]=0.
            v[j] = vitesse_faces(i_face, j);
 
          for (j = 0; j < dim; j++) // en 2D, xgr[3]=0.
            xgr[j] = centre_faces(i_face, j) - c[j];
 
          if (option_ != ENTRAINEMENT_SEUL)
            {
              src[0] += -2. * (w[1] * v[2] - w[2] * v[1]);
              src[1] += -2. * (w[2] * v[0] - w[0] * v[2]);
              src[2] += -2. * (w[0] * v[1] - w[1] * v[0]);
              src[0] += dw[2] * xgr[1] - dw[1] * xgr[2];
              src[1] += dw[0] * xgr[2] - dw[2] * xgr[0];
              src[2] += dw[1] * xgr[0] - dw[0] * xgr[1];
            }
          if (option_ != CORIOLIS_SEUL)
            {
              src[0] += w[2]*w[2]*xgr[0] - w[0]*w[2]*xgr[2] - w[0]*w[1]*xgr[1] + w[1]*w[1]*xgr[0];
              src[1] += w[0]*w[0]*xgr[1] - w[0]*w[1]*xgr[0] - w[1]*w[2]*xgr[2] + w[2]*w[2]*xgr[1];
              src[2] += w[1]*w[1]*xgr[2] - w[1]*w[2]*xgr[1] - w[2]*w[0]*xgr[0] + w[0]*w[0]*xgr[2];
            }
        }
      for (j = 0; j < dim; j++)
        acceleration_aux_faces(i_face, j) = src[j];
    }
  return acceleration_aux_faces;
}
 
/*! @brief resu=0; ajouter(resu); (calls ajouter() of the derived class)
 *
 */
DoubleTab& Terme_Source_Acceleration::calculer(DoubleTab& resu) const
{
  resu = 0;
  ajouter(resu);
  resu.echange_espace_virtuel();
  return resu;
}
 
/*! @brief Evaluates the acceleration and rotation fields at time t.
 *
 */
void Terme_Source_Acceleration::mettre_a_jour(double temps)
{
  //Cerr << "Terme_Source_Acceleration::mettre_a_jour temps=" << temps << finl;
  if (champ_acceleration_)
    champ_acceleration_->mettre_a_jour(temps);
  if (omega_)
    {
      omega_->mettre_a_jour(temps);
      domegadt_->mettre_a_jour(temps);
      centre_rotation_->mettre_a_jour(temps);
    }
  get_set_terme_source_post().mettre_a_jour(temps);
}
 
/*! @brief Overridden method from Source_base.
 *
 * Recognized keywords:
 *   "ACCELERATION" => returns terme_source_post_, homogeneous to d/dt(rho*v)
 *   Beware of the side effect of ajouter() (see comments of ajouter())
 *
 */
int Terme_Source_Acceleration::a_pour_Champ_Fonc(const Motcle& mot,
                                                 OBS_PTR(Champ_base) & ch_ref) const
{
  int ok = 0;
  if (mot == "ACCELERATION_TERME_SOURCE")
    {
      ch_ref = terme_source_post_.valeur();
      ok = 1;
    }
  return ok;
}
 
/*! @brief Returns eq_hydraulique_ !
 *
 */
const Navier_Stokes_std& Terme_Source_Acceleration::get_eq_hydraulique() const
{
  const Navier_Stokes_std& ns = ref_cast(Navier_Stokes_std, equation());
  return ns;
}