CDMATH/tests/examples/transport2d_ns/main.cxx

   1 //============================================================================
   2 // Author      : Anouar MEKKAS
   3 // Version     :
   4 // Description : Equation de transport lineaire 2D non structure
   5 //============================================================================
   6
   7 #include <iostream>
   8 #include <string>
   9 #include <cmath>
  10
  11 #include "Mesh.hxx"
  12 #include "Cell.hxx"
  13 #include "Face.hxx"
  14 #include "Field.hxx"
  15
  16 using namespace std;
  17
  18
  19 void conditions_initiales(Field& yField)
  20 {
  21     double rayon=0.15;
  22     double xcentre=0.25;
  23     double ycentre=0.25;
  24     Mesh myMesh=yField.getMesh();
  25     int nbCells=myMesh.getNumberOfCells();
  26     for (int j=0 ; j<nbCells ; j++)
  27     {
  28         double x = myMesh.getCell(j).x() ;
  29         double y = myMesh.getCell(j).y() ;
  30         double valX=(x-xcentre)*(x-xcentre);
  31         double valY=(y-ycentre)*(y-ycentre);
  32         double val=sqrt(valX+valY);
  33         if (val<rayon)
  34             yField(j) = 1.0;
  35         else
  36             yField(j) = 0.0;
  37     }
  38 }
  39
  40 void sigma_flux(double VitesseX, double VitesseY, double cfl, const Field& yField, const IntTab indexFacesPerio, double& dt, Field& SumFlux)
  41 {
  42     // Calculation of fluxes
  43     Mesh myMesh=yField.getMesh();
  44     int nbCells=myMesh.getNumberOfCells();
  45     double normU=sqrt(VitesseX*VitesseX+VitesseY*VitesseY);
  46     for (int j=0 ; j<nbCells ; j++)
  47     {
  48         Cell Cj=myMesh.getCell(j);
  49         int nbFace=Cj.getNumberOfFaces();
  50         double SumF=0.0;
  51         double minlengthFk=1.E30;
  52
  53         int cellCourante,cellAutre;
  54         for (int k=0 ; k<nbFace ; k++)
  55         {
  56             int indexFace=Cj.getFacesId()[k];
  57             Face Fk=myMesh.getFace(indexFace);
  58             double NormalX=Cj.getNormalVector(k,0);
  59             double NormalY=Cj.getNormalVector(k,1);
  60             double LengthFk = Fk.getMeasure();
  61             double UN=VitesseX*NormalX+VitesseY*NormalY;
  62
  63             minlengthFk=min(minlengthFk,LengthFk/fabs(UN));
  64             minlengthFk=min(minlengthFk,LengthFk/fabs(VitesseX));
  65             minlengthFk=min(minlengthFk,LengthFk/fabs(VitesseY));
  66
  67             double conc=0.0;
  68             cellCourante=j;
  69             cellAutre=-1;
  70
  71             if (!Fk.isBorder())
  72             {
  73                 int indexC1=Fk.getCellsId()[0];
  74                 int indexC2=Fk.getCellsId()[1];
  75                 /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
  76                 if ( indexC1 == j )
  77                 {
  78                     /* hypothese verifie */
  79                     cellCourante=indexC1;
  80                     cellAutre=indexC2;
  81                 } else if ( indexC2 == j )
  82                 {
  83                     /* hypothese non verifie */
  84                     cellCourante=indexC2;
  85                     cellAutre=indexC1;
  86                 }
  87                 // definir la cellule gauche et droite par le prduit vitesse * normale sortante
  88                 // si u*n>0 : rien a faire sinon inverser la gauche et la droite
  89                 if (UN>1.E-15)
  90                     conc=yField(cellCourante);
  91                 else
  92                     conc=yField(cellAutre);
  93             }else
  94             {
  95                 /* conditions aux limites neumann homogene */
  96                 if (Fk.getGroupName().compare("GAUCHE")==0 || Fk.getGroupName().compare("DROITE")==0)
  97                 {
  98                     if (UN>1.E-15)
  99                         conc=yField(cellCourante);
 100                     else
 101                         conc=0.0;
 102                 }
 103                 /* conditions aux limites periodiques */
 104                 if (Fk.getGroupName().compare("BAS")==0 || Fk.getGroupName().compare("HAUT")==0)
 105                 {
 106                         int indexFP=indexFacesPerio[indexFace];
 107                         /* une autre manière de recuperer l'index de la face periodique */
 108                         //int indexFP=M.getIndexFacePeriodic(indexFace);
 109                         Face Fp=myMesh.getFace(indexFP);
 110                         int indexCp=Fp.getCellsId()[0];
 111                         if (UN>1.E-15)
 112                             conc=yField(cellCourante);
 113                         else
 114                             conc=yField(indexCp);
 115                 }
 116             }
 117             SumF=SumF+UN*LengthFk*conc;
 118           }
 119         dt=cfl*minlengthFk/normU;
 120         SumFlux(j)=dt*SumF/Cj.getMeasure();
 121     }
 122 }
 123
 124 void EquationTransport2D(double tmax, double VitesseX, double VitesseY, double cfl, int freqSortie, const Mesh& myMesh, const string file)
 125 {
 126     /* Initial conditions */
 127     cout << "Construction of the initial condition …" << endl;
 128     Field yField("Y field",CELLS,myMesh,1) ;
 129     conditions_initiales(yField);
 130
 131     /*
 132      * MED output of the initial condition at t=0 and iter = 0
 133      */
 134     int iter=0;
 135     double time=0.;
 136     cout << "Saving the solution at T=" << time << "…" << endl;
 137     yField.setTime(time,iter);
 138     yField.writeMED(file);
 139     yField.writeVTK(file);
 140     yField.writeCSV(file);
 141     /* --------------------------------------------- */
 142
 143     /* Time loop */
 144     cout << "Resolution of the transport equation with an UPWIND scheme …" << endl;
 145     int ntmax=3;
 146     double dt;
 147     IntTab indexFacesPerio=myMesh.getIndexFacePeriodic();
 148     while (iter<ntmax && time <= tmax )
 149     {
 150         Field SumFlux("Fluxes sum",CELLS,myMesh,1) ;
 151         sigma_flux(VitesseX,VitesseY,cfl,yField,indexFacesPerio,dt,SumFlux);
 152         cout << "-- Iter: " << iter << ", Time: " << time << ", dt: " << dt << endl;
 153
 154         /* Advancing time step */
 155         yField-=SumFlux;
 156
 157         time+=dt;
 158         iter+=1;
 159         // Output every freq iterations
 160         if (iter%freqSortie==0)
 161         {
 162             yField.setTime(time,iter);
 163             yField.writeMED(file,false);
 164             yField.writeVTK(file,false);
 165             yField.writeCSV(file);
 166         }
 167     }
 168 }
 169
 170 int main()
 171 {
 172     cout << "Resolution of the 2D transport equation:" << endl;
 173     cout << "- DOMAIN: SQUARE" << endl;
 174     cout << "- MESH: TRIANGULAR, GENERATED WITH SALOME" << endl;
 175     cout << "- PERIODIC BC ON TOP AND BOTTOM" << endl;
 176     cout << "- HOMOGENEOUS NEUMANN BC ON LEFT AND RIGHT" << endl;
 177
 178     // Problem data
 179     double cfl=0.4;
 180     double VitesseX=1.0;
 181     double VitesseY=1.0;
 182     double tmax=1.;
 183     int freqSortie=10;
 184
 185     cout << "Construction of Cartesian mesh…" << endl;
 186     Mesh myMesh("../../tests/ressources/meshSquare.med");
 187     string fileOutPut="Exercice2";
 188     EquationTransport2D(tmax,VitesseX,VitesseY,cfl,freqSortie,myMesh,fileOutPut);
 189     cout << "CDMATH calculation done." << endl;
 190
 191     return 0;
 192 }