CDMATH/tests/examples/transport2d_s/main.cxx

   1 //============================================================================
   2 // Author      : Anouar MEKKAS
   3 // Version     :
   4 // Description : 2D linear transport equation on cartesian grid
   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 initial_conditions(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     /* Fluxes calculation */
  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                 /* Hypothesis: the cell of index indexC1 is the current cell index j */
  76                 if ( indexC1 == j )
  77                 {
  78                     /* Hypothesis verified */
  79                     cellCourante=indexC1;
  80                     cellAutre=indexC2;
  81                 } else if ( indexC2 == j )
  82                 {
  83                     /* Hypothesis not verified */
  84                     cellCourante=indexC2;
  85                     cellAutre=indexC1;
  86                 }
  87                 // Define left and right cells with the product velocity * outward normal
  88                 // If u*n>0, then nothing to do, else invert left and right
  89                 if (UN>1.E-15)
  90                     conc=yField(cellCourante);
  91                 else
  92                     conc=yField(cellAutre);
  93             }else
  94             {
  95                 /* Homogeneous Neumann boundary conditions */
  96                 if (Fk.getGroupName().compare("LeftEdge")==0 || Fk.getGroupName().compare("RightEdge")==0)
  97                 {
  98                     if (UN>1.E-15)
  99                         conc=yField(cellCourante);
 100                     else
 101                         conc=0.0;
 102                 }
 103                 /* Periodic boundary conditions */
 104                 if (Fk.getGroupName().compare("BottomEdge")==0 || Fk.getGroupName().compare("TopEdge")==0)
 105                 {
 106                         int indexFP=indexFacesPerio[indexFace];
 107                         /* une autre manière de recuperer l'index de la face periodique */
 108                         //int indexFP=myMesh.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     initial_conditions(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("Sum Flux",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         // Ouput 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 2D TRANSPORT EQUATION:" << endl;
 173     cout << "- DOMAIN: SQUARE" << endl;
 174     cout << "- MESH: CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << 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 a cartesian mesh …" << endl;
 186     double xinf=0.0;
 187     double xsup=1.0;
 188     double yinf=0.0;
 189     double ysup=1.0;
 190     int nx=100;
 191     int ny=100;
 192     Mesh myMesh(xinf,xsup,nx,yinf,ysup,ny);
 193     double eps=1.E-10;
 194     myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
 195     myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
 196     myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
 197     myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
 198     string fileOutPutCart="Exercice1";
 199     EquationTransport2D(tmax,VitesseX,VitesseY,cfl,freqSortie,myMesh,fileOutPutCart);
 200     cout << "CDMATH calculation done." << endl;
 201     return 0;
 202 }