Terme_Source_Qdm_lambdaup_VEF_Face.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 <Terme_Source_Qdm_lambdaup_VEF_Face.h>
#include <Domaine_VEF.h>
#include <Navier_Stokes_std.h>
#include <Probleme_base.h>
#include <Schema_Temps_base.h>
#include <Champ_P1NC.h>
#include <SFichier.h>
 
Implemente_instanciable(Terme_Source_Qdm_lambdaup_VEF_Face,"Source_Qdm_lambdaup_VEF_P1NC",Source_base);
// XD source_qdm_lambdaup source_base source_qdm_lambdaup BRACE This source term is a dissipative term which is intended
// XD_CONT to minimise the energy associated to non-conformscales u\' (responsible for spurious oscillations in some
// XD_CONT cases). The equation for these scales can be seen as: du\'/dt= -lambda. u\' + grad P\' where -lambda. u\'
// XD_CONT represents the dissipative term, with lambda = a/Delta t For Crank-Nicholson temporal scheme, recommended
// XD_CONT value for a is 2. NL2 Remark : This method requires to define a filtering operator.
 
Sortie& Terme_Source_Qdm_lambdaup_VEF_Face::printOn(Sortie& s ) const
{
  return s << que_suis_je() ;
}
 
Entree& Terme_Source_Qdm_lambdaup_VEF_Face::readOn(Entree& is )
{
  Motcle accouverte = "{" , accfermee = "}" ;
  Motcle motlu;
  is >> motlu;
  if(motlu!=accouverte)
    {
      Cerr << motlu << " --> " << finl;
      Cerr << "{ expected when reading lambda uprime source term " << finl;
      Cerr << "The syntax of keyword Source_Qdm_lambdaup has changed, see the" << finl;
      Cerr << "user manual of version 1.5.5 or later." << finl;
      Cerr << "The value of lambda, which can now be variable," << finl;
      Cerr << "must be preceded by a keyword:" << finl;
      Cerr << "Thus, Source_Qdm_lambdaup value becomes:" <<finl;
      Cerr << "Source_Qdm_lambdaup { lambda value }" <<finl;
      Process::exit();
    }
  Motcles les_mots(4);
  {
    les_mots[0] = "lambda";             // XD attr lambda floattant lambda_u REQ value of lambda
    les_mots[1] = "lambda_min";         // XD attr lambda_min floattant lambda_min OPT value of lambda_min
    les_mots[2] = "lambda_max";         // XD attr lambda_max floattant lambda_max OPT value of lambda_max
    les_mots[3] = "ubar_umprim_cible";  // XD attr ubar_umprim_cible floattant ubar_umprim_cible OPT value of
    // XD_CONT ubar_umprim_cible
  }
  lambda=-1;
  while(motlu!=accfermee)
    {
      is>>motlu;
      int rang=les_mots.search(motlu);
      switch(rang)
        {
        case 0:
          {
            is >> lambda;
            lambda_min=lambda;
            lambda_max=lambda;
            break;
          }
        case 1:
          {
            is >> lambda_min;
            break;
          }
        case 2:
          {
            is >> lambda_max;
            break;
          }
        case 3:
          {
            is >> cible;
            break;
          }
        default :
          {
            if(motlu!=accfermee)
              {
                Cerr << motlu
                     << " not recognised in lambda uprime source term "
                     << finl;
                Process::exit();
              }
          }
        }
    }
  if(lambda < 0)
    lambda=lambda_min;
  if (lambda < 0.)
    {
      Cerr << "Error reading lambda in " << que_suis_je() << finl;
      Cerr << " lambda must be defined and greater than 0 " << finl;
      exit();
    }
  Cerr << "Exiting readOn " << finl;
  return is ;
}
 
void Terme_Source_Qdm_lambdaup_VEF_Face::associer_pb(const Probleme_base& pb)
{
  int nb_eqn = pb.nombre_d_equations();
  int ok=0;
  for (int i=0; i<nb_eqn; i++)
    {
      const Equation_base& eqn = pb.equation(i);
      if  (sub_type(Navier_Stokes_std,eqn))
        {
          la_vitesse = ref_cast(Champ_P1NC,eqn.inconnue());
          associer_domaines(eqn.domaine_dis(),eqn.domaine_Cl_dis());
          i = nb_eqn;
          ok = 1;
        }
      else
        Cerr << "Lambdaup : " << eqn.que_suis_je() << finl;
    }
 
  if (!ok)
    {
      Cerr << "TRUST error in " << que_suis_je()  << finl;
      Cerr << "No hydraulic equation found in the problem" << finl;
      exit();
    }
}
 
void Terme_Source_Qdm_lambdaup_VEF_Face::associer_domaines(const Domaine_dis_base& domaine_dis,
                                                           const Domaine_Cl_dis_base& domaine_Cl_dis)
{
  le_dom_VEF = ref_cast(Domaine_VEF, domaine_dis);
  le_dom_Cl_VEF = ref_cast(Domaine_Cl_VEF, domaine_Cl_dis);
}
 
 
DoubleTab& Terme_Source_Qdm_lambdaup_VEF_Face::ajouter(DoubleTab& resu) const
{
  static double rapport_old=1.;
  const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
  const DoubleVect& volumes_entrelaces=domaine_VEF.volumes_entrelaces();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
  const DoubleTab& vitesse=la_vitesse->valeurs();
  DoubleTab ubar(vitesse);
  DoubleTab uprime(vitesse);
  const int nb_faces = domaine_VEF.nb_faces();
  int face, k;
 
  int nb_comp = resu.line_size();
 
  la_vitesse->filtrer_L2(ubar);
  uprime-=ubar;
  // Optimization: combine 2 mp_norme_vect into 1 collective call
  double normbar_carre = local_carre_norme_vect(ubar);
  double normprim_carre = local_carre_norme_vect(uprime);
  Process::mp_sum_for_each(normbar_carre, normprim_carre);
  double normbar = sqrt(normbar_carre);
  double normprim = sqrt(normprim_carre);
  //double rapport=normbar/(normprim+DMINFLOAT) ;
  double rapport=(normprim!=0?normbar/normprim:0);
  double& l=lambda;
  if( (rapport < cible) && (rapport<rapport_old) ) l*=2.;
  if( (rapport < cible) && (rapport>rapport_old) ) l*=1.1;
  if( (rapport>cible) && (rapport>rapport_old)  )  l*=0.95;
  l=std::min(lambda, lambda_max);
  l=std::max(lambda_min, lambda);
  Cerr << "|ubar| = " << normbar << finl;
  Cerr << "|ubar|/|uprim| = " << rapport << finl;
  Cerr << "lambda " << lambda << finl;
  for(face=0; face<nb_faces; face++)
    for(k=0; k< dimension; k++)
      {
        double v=volumes_entrelaces(face)*porosite_face(face);
        if(nb_comp>1)
          {
            uprime(face,k)*=v;
          }
        else
          {
            uprime(face)*=v;
          }
      }
 
  double dt=la_vitesse->equation().schema_temps().pas_de_temps();
  resu.ajoute(-lambda/dt, uprime);
  rapport_old=rapport;
  return resu;
}
 
DoubleTab& Terme_Source_Qdm_lambdaup_VEF_Face::calculer(DoubleTab& resu) const
{
  resu = 0;
  return ajouter(resu);
}
 
void Terme_Source_Qdm_lambdaup_VEF_Face::mettre_a_jour(double temps)
{
  ;
}
 
double Terme_Source_Qdm_lambdaup_VEF_Face::norme_H1(const DoubleTab& vitesse) const
{
  const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
  const Domaine& domaine = domaine_VEF.domaine();
 
  double dnorme_H1,norme_H1_comp,int_grad_elem,norme_grad_elem;
  int face_globale;
 
  //Compute the H1 norm of a P1NC unknown.
  //The algorithm accounts for the following constraints:
  //- the unknown may have multiple components
  //  (i.e. be a scalar or a vector)
  //- the problem dimension is arbitrary (1, 2 or 3).
  //WARNING: prisms are not supported.
 
  dnorme_H1=0.;
  for (int composante=0; composante<vitesse.line_size(); composante++)
    {
      norme_H1_comp=0.; //to avoid accumulation
      for (int K=0; K<domaine.nb_elem(); K++) //loop over elements
        {
          norme_grad_elem=0.; //to avoid accumulation
          for (int i=0; i<dimension; i++) //loop over problem dimension
            {
              int_grad_elem=0.; //to avoid accumulation
              for (int face=0; face<domaine.nb_faces_elem(); face++) //loop over faces of element "K"
                {
                  face_globale = domaine_VEF.elem_faces(K,face);
 
                  int_grad_elem += vitesse(face_globale,composante)*
                                   domaine_VEF.face_normales(face_globale,i)*
                                   domaine_VEF.oriente_normale(face_globale,K);
                } //end of for loop over "face"
 
              norme_grad_elem += int_grad_elem*int_grad_elem;
            } //end of for loop over "i"
 
          norme_H1_comp += norme_grad_elem/domaine_VEF.volumes(K);
        } //end of for loop over "K"
 
      dnorme_H1 += norme_H1_comp;
    } // end of for loop over "composante"
 
  return sqrt(dnorme_H1);
}
 
double Terme_Source_Qdm_lambdaup_VEF_Face::norme_L2_H1(const DoubleTab& u) const
{
  double pas_de_temps = la_vitesse->equation().schema_temps().pas_de_temps();
 
  return ((1./pas_de_temps)*norme_L2(u))+norme_H1(u);
}
 
double Terme_Source_Qdm_lambdaup_VEF_Face::norme_L2(const DoubleTab& u) const
{
  const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
  const DoubleVect& volumes = domaine_VEF.volumes_entrelaces();
  const int nb_faces = domaine_VEF.nb_faces();
 
  int i=0;
  double norme =0;
  int nb_compo_=u.line_size();
 
  for(; i<nb_faces; i++)
    {
      double s=0;
      for(int j=0; j<nb_compo_; j++)
        s+=u(i,j)*u(i,j);
      norme+=volumes(i)*s;
    }
 
  return sqrt(norme);
}