Ch_front_Vortex.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 <Ch_front_Vortex.h>
#include <Domaine.h>
#include <LecFicDiffuse.h>
#include <Interprete.h>
#include <Domaine_VF.h>
#include <time.h>
#include <EcrFicCollecte.h>
 
Implemente_instanciable_sans_destructeur(Ch_front_Vortex,"Champ_front_Vortex",Champ_front_var_instationnaire);
 
 
Ch_front_Vortex::~Ch_front_Vortex()
{
  sauvegarder_vortex();
}
 
/*! @brief Prints the field to an output stream.
 *
 * Prints the size of the field and the (constant) value on
 *     the boundary.
 *
 * @param (Sortie& os) an output stream
 * @return (Sortie&) the modified output stream
 */
Sortie& Ch_front_Vortex::printOn(Sortie& os) const
{
  const DoubleTab& tab=valeurs();
  os << tab.size() << " ";
  for(int i=0; i<tab.size(); i++)
    os << tab(0,i);
  return os;
}
 
/*! @brief Reads the field from an input stream.
 *
 * Format:
 *       Champ_front_Vortex nb_compo vrel_1 ... [vrel_i]
 *
 * @param (Entree& is) an input stream
 * @return (Entree& is) the modified input stream
 */
Entree& Ch_front_Vortex::readOn(Entree& is)
{
  if(dimension!=3)
    {
      Cerr << " Vortex Boundary condition requests 3D problem " << finl;
      exit();
    }
 
  fixer_nb_comp(dimension);
 
  Nom nom;
  is >> nom >> geom >> nu >> utau;
 
  mon_domaine = ref_cast(Domaine, Interprete::objet(nom));
 
  return is;
}
 
/*! @brief Returns the object upcast to Champ_front_base&
 *
 * @param (Champ_front_base& ch)
 * @return (Champ_front_base&) (*this) upcast to Champ_front_base&
 */
Champ_front_base& Ch_front_Vortex::affecter_(const Champ_front_base& ch)
{
  return *this;
}
 
 
int Ch_front_Vortex::my_rand( )
{
  if (first_rand == 1)
    {
      srand((int)time (nullptr));
      first_rand = 0;
    }
  return (rand ());
}
 
void Ch_front_Vortex::sauvegarder_vortex()
{
  Nom fichier = "vortex.sauv";
  if(la_frontiere_dis)
    {
      Cerr << "Saving vortices in " << fichier << finl;
      EcrFicCollecte fic(fichier);
      fic.setf(ios::scientific);
 
      int nbvortex =  xvort.size();
 
      fic << temps << finl;
      fic << nbvortex << finl;
 
      for (int i=0; i<nbvortex; i++)
        {
          fic << xvort(i) << " " << yvort(i) << " " << zvort(i) << " " << tvort(i) << " " ;
          fic << svort(i) << " " << fvort[i] << " " << gamma(i) << " " << sigma(i) << finl;
        }
      fic.close();
    }
}
 
void Ch_front_Vortex::reprendre_vortex()
{
  Cerr << " Ch_front_Vortex::reprendre_vortex" << finl;
 
  Nom fichier = "vortex.sauv";
  if (nproc()>1) fichier=fichier.nom_me(me());
  if(la_frontiere_dis)
    {
      Cerr << "Resumption of vortices in " << fichier << finl;
      EFichier fic(fichier);
      fic.setf(ios::scientific);
 
      if(fic.fail())
        {
          Cerr << "Unable to open file " << fichier << finl;
          exit();
        }
 
      double tps;
      int nbvortex;
 
      fic >> tps;
      fic >> nbvortex;
 
      if(!(std::fabs(temps-tps)/(std::fabs(tps)+1.e-24)<1.e-6))
        {
          Cerr << " time : " << tps << " in " << fichier << " does not match the time of resumed " << temps << finl;
          exit();
        }
 
      Cerr << " resumed " << nbvortex << " vortex at time " << temps << finl;
 
      for (int i=0; i<nbvortex; i++)
        {
          fic >> xvort(i) >> yvort(i) >> zvort(i) >> tvort(i);
          fic >> svort(i) >> fvort[i] >> gamma(i) >> sigma(i);
        }
 
      Cerr << " resumed made " << finl;
    }
}
 
 
int Ch_front_Vortex::initialiser(double un_temps, const Champ_Inc_base& inco)
{
  if (!Champ_front_var_instationnaire::initialiser(un_temps,inco))
    return 0;
 
  Domaine& domaine=mon_domaine.valeur();
  const Frontiere_dis_base& fr_dis=frontiere_dis();
  const Frontiere& frontiere=fr_dis.frontiere();
  const int nb_faces=frontiere.nb_faces();
  const int ndeb=frontiere.num_premiere_face();
  const int nfin=ndeb+nb_faces;
  const Faces& faces=frontiere.faces();
  int nbsf=faces.nb_som_faces();
  const Domaine_VF& zvf = ref_cast(Domaine_VF, domaine_dis());
  const DoubleTab& xv=zvf.xv();
 
  double x,y,z;
  double dx,dy,dz;
  double norme;
 
  double  C_mu = 0.09;
 
 
  u.resize(nb_faces);
  v.resize(nb_faces);
  w.resize(nb_faces);
  u_moy.resize(nb_faces);
  dudy.resize(nb_faces);
  k.resize(nb_faces);
  eps.resize(nb_faces);
  Wk.resize(nb_faces);
 
  for(int face=ndeb; face<nfin; face++)
    {
      Wk(face) = 2.*drand48()-1.;
    }
 
 
  /////////////////////////////////////////////////////////
  // geometric characteristics associated with the boundary
  /////////////////////////////////////////////////////////
 
  if(geom=="circle")
    {
 
      double x_xmin=1.e6;
      double x_xmax=-1.e6;
      double y_xmin=1.e6;
      double y_xmax=-1.e6;
      double z_xmin=1.e6;
      double z_xmax=-1.e6;
      double x_ymin=1.e6;
      double x_ymax=-1.e6;
      double y_ymin=1.e6;
      double y_ymax=-1.e6;
      double z_ymin=1.e6;
      double z_ymax=-1.e6;
 
 
      for(int face=ndeb; face<nfin; face++)
        {
          for(int kk=0; kk<nbsf; kk++)
            {
              x=domaine.coord(faces.sommet(face,kk),0);
              y=domaine.coord(faces.sommet(face,kk),1);
              z=domaine.coord(faces.sommet(face,kk),2);
 
              if(y<y_ymin)
                {
                  y_ymin=y;
                  x_ymin=x;
                  z_ymin=z;
                }
              if(y>y_ymax)
                {
                  y_ymax=y;
                  x_ymax=x;
                  z_ymax=z;
                }
 
              if(x<x_xmin)
                {
                  x_xmin=x;
                  y_xmin=y;
                  z_xmin=z;
                }
              if(x>x_xmax)
                {
                  x_xmax=x;
                  y_xmax=y;
                  z_xmax=z;
                }
            }
        }
 
      if(!est_egal(y_ymin,y_ymax))
        {
          dx=x_ymax-x_ymin;
          dy=y_ymax-y_ymin;
          dz=z_ymax-z_ymin;
 
          Ox=0.5*(x_ymax+x_ymin);
          Oy=0.5*(y_ymax+y_ymin);
          Oz=0.5*(z_ymax+z_ymin);
        }
      else
        {
          dx=x_xmax-x_xmin;
          dy=y_xmax-y_xmin;
          dz=z_xmax-z_xmin;
 
          Ox=0.5*(x_xmax+x_xmin);
          Oy=0.5*(y_xmax+y_xmin);
          Oz=0.5*(z_xmax+z_xmin);
        }
 
      R=0.5*sqrt(dx*dx+dy*dy+dz*dz);
 
      surf = M_PI*R*R;
 
      Cerr << " Coordinates associated with the boundary : " << Ox << " " << Oy << " " << Oz << finl;
      Cerr << " Radius computed : " << R << finl;
 
    }
  else if(geom=="channel")
    {
      //
      // the channel is assumed to be oriented along x,
      // periodic along z (and x)
      // the inlet boundary centered at the origin
      // with height and width: 2
      // i.e. (x,y,z) : [0;0]x[-1;+1]x[-1;+1]
 
      R=1.;
 
      Ox=0.;
      Oy=0.;
      Oz=0.;
 
      surf = 4.;
    }
 
  ///////////////////////////////////////////////////////////////
  // local frame (t1;t2) associated with the boundary (assumed to be planar)
  ///////////////////////////////////////////////////////////////
  if (nb_faces==0)
    {
      nb_vortex=0;
    }
  else
    {
      nx=-zvf.face_normales(ndeb,0);  // normals oriented towards the interior of the domain
      ny=-zvf.face_normales(ndeb,1);
      nz=-zvf.face_normales(ndeb,2);
 
      norme = sqrt( nx*nx + ny*ny + nz*nz );
 
      nx/=norme;
      ny/=norme;
      nz/=norme;
 
      t1x = -ny;
      t1y =  nx;
      t1z =  0.;
 
      if( (est_egal(nx,0.)) && (est_egal(ny,0.)) )
        {
          t1x =  0.;
          t1y =  nz;
          t1z = -ny;
        }
 
      norme = sqrt( t1x*t1x + t1y*t1y + t1z*t1z );
 
      t1x/=norme;
      t1y/=norme;
      t1z/=norme;
 
      t2x =  ny*t1z - nz*t1y;
      t2y =  nz*t1x - nx*t1z;
      t2z =  nx*t1y - ny*t1x;
      /*
           Cerr << " t1x " << t1x << finl;
           Cerr << " t1y " << t1y << finl;
           Cerr << " t1z " << t1z << finl;
           Cerr << " t2x " << t2x << finl;
           Cerr << " t2y " << t2y << finl;
           Cerr << " t2z " << t2z << finl;
      */
      ///////////////////////////////////////////////////////////////////
      // initialization of the fields: velocity - k - epsilon (analytical)
      ///////////////////////////////////////////////////////////////////
 
      double d,dplus,Reytau;
      double kappa = 0.41;
      double beta = 5.1;
      double tmp=15*15*15*15;
      for(int face=ndeb; face<nfin; face++)
        {
          x=xv(face,0);
          y=xv(face,1);
          z=xv(face,2);
 
          dx=x-Ox;
          dy=y-Oy;
          dz=z-Oz;
 
          if(geom=="circle")
            {
              d = R - sqrt(dx*dx+dy*dy+dz*dz) ;  // distance to the wall
            }
          else
            {
              d = R - std::fabs(y) ;  // distance to the wall
            }
 
          if(d<=0.)
            {
              Cerr << " the computation of a zero or negative distance is prejudicial for the rest of the calculation" << finl;
              Cerr << " Problem certainly related to the calculation of the size of domain or origin related to the bondary" << finl;
              exit();
            }
 
          dplus   = d*utau/nu;
          Reytau  = R*utau/nu;
 
          if(dplus<11.)
            {
              u_moy(face) = utau*dplus ;
            }
          else
            {
              u_moy(face) = utau*((1./kappa)*log(dplus)+beta) ;
            }
 
          if(dplus<11.)
            {
              dudy(face) = utau*utau/nu ;  // -> discontinuous dudy gradient at dplus = 11
            }
          else
            {
              dudy(face) = utau/(kappa*d) ;  // -> generates a discontinuity on u' at this location
            }
 
          k(face)   = utau*utau * ( 0.07*dplus*dplus*exp(-dplus/8)+4.5*(1.-exp(-dplus/20.))/(1.+4.*dplus/Reytau) ) ;
          eps(face) = (utau*utau*utau*utau/nu) * (1./(0.41 *(pow(dplus*dplus*dplus*dplus+tmp,0.25)))) ;
        }
 
      u = u_moy;
 
      //////////////////////////////////
      // evaluation of the number of vortices
      /////////////////////////////////
 
      double sigma_moy = 0.;
 
      for(int face=ndeb; face<nfin; face++)
        {
          double sx=zvf.face_normales(face,0);
          double sy=zvf.face_normales(face,1);
          double sz=zvf.face_normales(face,2);
 
          sigma_moy += sqrt(sx*sx+sy*sy+sz*sz) * pow(k(face),1.5) / eps(face);
 
        }
      sigma_moy *= pow(C_mu,0.75);
      sigma_moy /= surf;
      // We cannot use modf (see Schema_Temps_base::limpr) because
      // nb_vortex is used to size arrays!
      nb_vortex = (int) ( surf / (M_PI * sigma_moy *sigma_moy) );
    }
  Cerr << "nb_vortex " <<nb_vortex << finl;
 
  //nb_vortex=1;
  //nb_vortex=nb_faces;
 
  xvort.resize(nb_vortex);
  yvort.resize(nb_vortex);
  zvort.resize(nb_vortex);
  tvort.resize(nb_vortex);
  svort.resize(nb_vortex);
  fvort.resize_array(nb_vortex);
  gamma.resize(nb_vortex);
  sigma.resize(nb_vortex);
 
  //////////////////////////////////////////////////////
  // Computation of the minimum lifetime of the vortices
  //////////////////////////////////////////////////////
 
  double dt_vortex,dt_min=1.e6;
 
  for (int i=0; i<nb_faces; i++) // estimation of the required dt_min
    {
      dt_vortex=5.*C_mu*pow(k(i),3./2.)/(eps(i)*u_moy(i));
      if(dt_vortex<dt_min) dt_min=dt_vortex;
    }
 
  Cerr << " DT_MIN VORTEX " << dt_min << finl;
 
  /////////////////////////////////////////////////
  // vortex initialisation by random drawing
  /////////////////////////////////////////////////
 
  /*
           for (int i=0;i<nb_vortex;i++)  // different random at each new draw
           {
           int face = my_rand() % nb_faces;
           double rand1 = (my_rand() % 100)/100.;
           double rand2 = (my_rand() % 100)/100.;
           xvort(i)=xv(face,0);
           yvort(i)=xv(face,1);
           zvort(i)=xv(face,2);
           tvort(i)=rand1*5.*C_mu*pow(k(face),3./2.)/(eps(face)*u_moy(face));
           if (rand2<0.5) svort(i)=-1.;
           else svort(i)=1.;
 
           Cerr << " Vortex " << i << finl;
           Cerr << " Face " << face << " " << xvort(i) << " " <<yvort(i) << " " <<zvort(i) << finl;
           Cerr <<" " << tvort(i) << " " <<svort(i) << finl;
           Cerr << " " << finl;
 
           }
  */
 
  for (int i=0; i<nb_vortex; i++) // identical random sequence at each new draw
    {
      int face = (int)(nb_faces*drand48());
      xvort(i)=xv(face,0);
      yvort(i)=xv(face,1);
      zvort(i)=xv(face,2);
      tvort(i)=drand48()*5.*C_mu*pow(k(face),3./2.)/(eps(face)*u_moy(face));
      if (drand48()<0.5) svort(i)=-1.;
      else svort(i)=1.;
      /*
             Cerr << " Vortex " << i << finl;
             Cerr << " Face " << face << " " << xvort(i) << " " <<yvort(i) << " " <<zvort(i) << finl;
             Cerr <<" " << tvort(i) << " " <<svort(i) << finl;
             Cerr << " " << finl;
      */
    }
 
  mettre_a_jour(un_temps);
 
  return 1;
}
 
void  Ch_front_Vortex::deplacement_vortex(double dx,double dy,double dz,double asigma,double agamma,
                                          double& vv,double& ww)
{
  double norme2 = dx*dx + dy*dy + dz*dz ;
 
  if(norme2<(4.*asigma) && norme2>1.e-8)
    {
      double e1 = exp(-norme2/(2.*asigma*asigma));
      double f  = (1.-e1) * agamma * e1 / (2.* M_PI * norme2 );
 
      double dt1 = dx*t1x + dy*t1y + dz*t1z ;
      double dt2 = dx*t2x + dy*t2y + dz*t2z ;
 
      vv += dt2*f;
      ww -= dt1*f;
    }
}
 
void Ch_front_Vortex::mettre_a_jour(double tps)
{
  const Frontiere_dis_base& fr_dis=frontiere_dis();
  const Frontiere& frontiere=fr_dis.frontiere();
  const int nb_faces=frontiere.nb_faces();
  const int ndeb=frontiere.num_premiere_face();
  const int nfin=ndeb+nb_faces;
  DoubleTab& les_val=valeurs();
 
  const Domaine_VF& zvf = ref_cast(Domaine_VF, domaine_dis());
  const DoubleTab& xv=zvf.xv();
 
  double alpha=0;
  double dist,dist2;
  double x,yter,z;
  double xp,yp,zp;
  double xvort_mir,yvort_mir,zvort_mir;
  double xvort_per,yvort_per,zvort_per;
  double dx,dy,dz;
 
  double  C_mu = 0.09;
 
  if(init)
    {
      init=0;
      temps=tps;
 
      if(temps!=0.) // this is a restart
        reprendre_vortex();
    }
 
  //////////////////////////////////////////////////
  // vortex characteristics (position-vorticity)
  //////////////////////////////////////////////////
  if (nb_vortex!=0)
    alpha = 4.*sqrt(M_PI*surf/(3.*nb_vortex*(2.*log(3.)-3.*log(2.))));
 
  for (int i=0; i<nb_vortex; i++)
    {
      x=xvort(i);
      yter=yvort(i);
      z=zvort(i);
 
      dist=1.e6;
      int face=-1;
 
      for(int j=ndeb; j<nfin; j++)
        {
          dist2=(xv(j,0)-x)*(xv(j,0)-x)+(xv(j,1)-yter)*(xv(j,1)-yter)+(xv(j,2)-z)*(xv(j,2)-z);
 
          if(dist2<dist)
            {
              dist=dist2;
              face=j;
            }
        }
      gamma(i)     = alpha * sqrt(k(face)) * svort(i);
      sigma(i)     = pow(C_mu,0.75) * pow(k(face),1.5) / eps(face);
      fvort[i]     = face ;
    }
 
  ///////////////////////////////////////////////////////////////////////////////////
  // computation of the velocity induced by the vortices at each face of the boundary
  ///////////////////////////////////////////////////////////////////////////////////
 
  for(int face=ndeb; face<nfin; face++)
    {
      x=xv(face,0);
      yter=xv(face,1);
      z=xv(face,2);
 
      x -= Ox ;  // NOTE: For practical reasons,
      yter -= Oy ;  // point coordinates are expressed
      z -= Oz ;  // relative to the boundary origin Ox,Oy,Oz
 
      v(face) = 0.;  // velocity induced by vortices at the face along t1
      w(face) = 0.;  // velocity induced by vortices at the face along t2
 
      for (int i=0; i<nb_vortex; i++)
        {
          dx = (xvort(i)-Ox) - x;
          dy = (yvort(i)-Oy) - yter;
          dz = (zvort(i)-Oz) - z;
 
          deplacement_vortex(dx,dy,dz,sigma(i),gamma(i),v(face),w(face));
 
 
          ////////////////////////////////////
          // boundary condition treatment
          ////////////////////////////////////
 
          //***************************************
          // Walls: mirror vortex treatment
          //***************************************
 
          if(geom=="circle")
            {
              xp = x * R / (sqrt(x*x+yter*yter+z*z)) ;
              yp = yter * R / (sqrt(x*x+yter*yter+z*z)) ;
              zp = z * R / (sqrt(x*x+yter*yter+z*z)) ;
            }
          else //channel
            {
              xp = x;
              yp = R;
              if(yter<0) yp = -R;
              zp = z;
            }
 
          xvort_mir = 2. * xp - (xvort(i)-Ox);
          yvort_mir = 2. * yp - (yvort(i)-Oy);
          zvort_mir = 2. * zp - (zvort(i)-Oz);
 
          dx = xvort_mir - x;
          dy = yvort_mir - yter;
          dz = zvort_mir - z;
 
          deplacement_vortex(dx,dy,dz,sigma(i),gamma(i),v(face),w(face));
 
 
          //***************************************
          // Periodicity: planar channel only
          //***************************************
 
          if(geom=="channel")
            {
 
              xvort_per = (xvort(i)-Ox);
              yvort_per = (yvort(i)-Oy);
              if( (zvort(i)-Oz) >0.) zvort_per = -2.*R + (zvort(i)-Oz) ;
              else                            zvort_per =  2.*R + (zvort(i)-Oz);
 
              dx = xvort_per - x;
              dy = yvort_per - yter;
              dz = zvort_per - z;
 
              deplacement_vortex(dx,dy,dz,sigma(i),gamma(i),v(face),w(face));
 
            }
 
        }//loop over vortices
 
 
      ////////////////////////////////////////////////////////////////////////////
      // computation of the longitudinal velocity fluctuation via the Langevin equation
      ////////////////////////////////////////////////////////////////////////////
 
 
      //     double up = u(face) - u_moy(face);
      double up= u(face);
      if(tps==0.) up=0.;
      double vp;
 
      if(geom=="channel") vp = w(face);
      else //(geom=="circle")
        {
          double norme = sqrt(x*x+yter*yter+z*z);
          double costheta = ( x*t1x + yter*t1y + z*t1z ) / norme;
          double sintheta = ( x*t2x + yter*t2y + z*t2z ) / norme;
          vp = costheta*v(face) + sintheta*w(face) ;
        }
 
      ///// Wiener process: "Gaussian" random draw satisfying <dW>=0 and <dW^2>=1 /////////
      /////
      ///// http://www.taygeta.com/random/gaussian.html
 
      double x1,x2,ybis,Wkp1,dW;
 
      do
        {
          x1 = 2. * drand48() - 1. ;
          x2 = 2. * drand48() - 1. ;
          ybis = x1 * x1 + x2 * x2;
        }
      while ( ybis >= 1. );
 
      Wkp1     =  x1 * sqrt( (-1.0 * log( ybis ) ) / ybis );  // Working on dW -> to ensure <dW^2>=1, the coefficient must be modified
      dW       =  Wkp1 - Wk(face) ;                       // compared to the formulation proposed in 'taygeta'
      Wk(face) =  Wkp1 ;
 
      /*
        double Wkp1,dW;
 
        Wkp1     =  2.*drand48()-1.;
        dW       = (Wkp1 - Wk(face))*sqrt(3./2.) ;
        Wk(face) =  Wkp1 ;
 
 
        double W,Wp1,dW;
        W     =  2.*drand48()-1.;
        Wp1   =  2.*drand48()-1.;
        dW    = (Wp1 - W)*sqrt(3./2.) ;
      */
 
 
      //////////////////////////////////
 
      double C0,C1,C2;
 
      C0=14./15.;
      C1=1.8;
      C2=0.6;
 
      double T=k(face)/eps(face);
 
      double dt = tps - temps;
 
      if(dt!=0.) // implementation following Sergent's thesis + EDF paper (up to a factor 2/3 relative to C2)
        {
          up += dt*( -(C1/(2.*T))*up  + (C2-1.)*dudy(face)*vp + sqrt(C0*eps(face))*dW/dt );
 
          //        up += dt*( -(C1/(2.*T))*up  + ((2./3.)*C2-1.)*dudy(face)*vp + sqrt(C0*eps(face))*dW/dt );
        }
      /*
             if(dt!=0.) // implementation following what is coded in Saturne
             {
             up += dt*(-(1.-(2./3.)*C2)*dudy(face)*vp)  + 2.*sqrt((2./3.)*(C1-1.)*eps(face)*dt)*dW;
             up /= (1.+5.*C1*dt/T);
             }
      */
 
      //     u(face) = u_moy(face) + up;
      u(face) = up;
 
 
      les_val(face,0) = u(face)*nx + v(face)*t1x + w(face)*t2x ;
      les_val(face,1) = u(face)*ny + v(face)*t1y + w(face)*t2y ;
      les_val(face,2) = u(face)*nz + v(face)*t1z + w(face)*t2z ;
 
      /*
        les_val(face,0) =  gamma(face) ;
        les_val(face,1) =  sigma(face) ;
        les_val(face,2) =  fvort(face) ;
 
        les_val(face,0) =  k(face)  ;
        les_val(face,1) =  eps(face)  ;
        les_val(face,2) =  u_moy(face)  ;
 
      */
    }//loop over faces
 
 
  //////////////////////////
  // vortex displacement
  //////////////////////////
 
 
  double dt = tps - temps;
  temps=tps;
 
 
  for (int i=0; i<nb_vortex; i++)
    {
      int face = fvort[i];
 
      x = xvort(i) + dt * ( v(face)*t1x + w(face)*t2x ); // Vortex displacements only occur
      yter = yvort(i) + dt * ( v(face)*t1y + w(face)*t2y ); // in the inlet plane, hence the u.n
      z = zvort(i) + dt * ( v(face)*t1z + w(face)*t2z ); // components are excluded
 
      x -= Ox ;  // NOTE: For practical reasons,
      yter -= Oy ;  // point coordinates are expressed
      z -= Oz ;  // relative to the boundary origin Ox,Oy,Oz
 
 
      if( (geom=="circle") && ((x*x+yter*yter+z*z)>=(R*R)) ) //wall
        {
          // wall: do nothing - leave the vortex at its current position
        }
      else if( (geom=="channel") && (std::fabs(yter)>R) ) //wall
        {
          // wall: do nothing - leave the vortex at its current position
        }
      else if( (geom=="channel") && (z>R) ) //periodicity
        {
          xvort(i) = x + Ox;
          yvort(i) = yter + Oy;
          zvort(i) = -2.*R + z + Oz;
        }
      else if( (geom=="channel") && (z<-R) ) //periodicity
        {
          xvort(i) = x + Ox;
          yvort(i) = yter + Oy;
          zvort(i) = 2.*R - z + Oz;
        }
      else
        {
          xvort(i) = x + Ox;
          yvort(i) = yter + Oy;
          zvort(i) = z + Oz;
        }
    }
 
  ///////////////////////////////////////////////////////
  // end of vortex lifetime - generation of new vortex
  ///////////////////////////////////////////////////////
 
  for (int i=0; i<nb_vortex; i++)
    {
      tvort(i) -= dt ;
 
      if(tvort(i)<=0.)// end of vortex lifetime
        {
          /*
            int face = my_rand() % nb_faces;            // different random at each new draw
            //       double rand1 = (my_rand() % 100)/100.;
            double rand2 = (my_rand() % 100)/100.;
            xvort(i)=xv(face,0);
            yvort(i)=xv(face,1);
            zvort(i)=xv(face,2);
            tvort(i)=5.*C_mu*pow(k(face),3./2.)/(eps(face)*u_moy(face));
            if (rand2<0.5) svort(i)=-1.;
            else svort(i)=1.;
          */
          int face = (int)(nb_faces*drand48());       // identical random sequence at each execution
          xvort(i)=xv(face,0);
          yvort(i)=xv(face,1);
          zvort(i)=xv(face,2);
          tvort(i)=5.*C_mu*pow(k(face),3./2.)/(eps(face)*u_moy(face));
          if (drand48()<0.5) svort(i)=-1.;
          else svort(i)=1.;
 
          Cerr << " Time " << temps << " Vortex " << i << finl;
          Cerr << " Face " << face << " " << xvort(i) << " " <<yvort(i) << " " <<zvort(i) << finl;
          Cerr <<" " << tvort(i) << " " <<svort(i) << finl;
          Cerr << " " << finl;
        }
    }
}