CDMATH/tests/examples/transport2d_ns/main2.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 struct parameters {
  19     double cfl;
  20     double VitesseX;
  21     double VitesseY;
  22     int freqSortie;
  23     std::string outFile;
  24 };
  25
  26 void conditions_initiales(Field& yField)
  27 {
  28     double rayon=0.15;
  29     double xcentre=0.25;
  30     double ycentre=0.25;
  31     Mesh myMesh=yField.getMesh();
  32     int nbCells=myMesh.getNumberOfCells();
  33     for (int j=0 ; j<nbCells ; j++)
  34     {
  35         double x = myMesh.getCell(j).x() ;
  36         double y = myMesh.getCell(j).y() ;
  37         double valX=(x-xcentre)*(x-xcentre);
  38         double valY=(y-ycentre)*(y-ycentre);
  39         double val=sqrt(valX+valY);
  40         if (val<rayon)
  41             yField(j) = 1.0;
  42         else
  43             yField(j) = 0.0;
  44     }
  45 }
  46
  47 void sigma_flux(const parameters & p, const Field& yField, const IntTab indexFacesPerio, double& dt, Field& SumFlux)
  48 {
  49     /* Calcul des flux */
  50     Mesh myMesh=yField.getMesh();
  51     int nbCells=myMesh.getNumberOfCells();
  52     double normU=sqrt(p.VitesseX*p.VitesseX+p.VitesseY*p.VitesseY);
  53     for (int j=0 ; j<nbCells ; j++)
  54     {
  55         Cell Cj=myMesh.getCell(j);
  56         int nbFace=Cj.getNumberOfFaces();
  57         double SumF=0.0;
  58         double minlengthFk=1.E30;
  59
  60         int cellCourante,cellAutre;
  61         for (int k=0 ; k<nbFace ; k++)
  62         {
  63             int indexFace=Cj.getFacesId()[k];
  64             Face Fk=myMesh.getFace(indexFace);
  65             double NormalX=Cj.getNormalVector(k,0);
  66             double NormalY=Cj.getNormalVector(k,1);
  67             double LengthFk = Fk.getMeasure();
  68             double UN=p.VitesseX*NormalX+p.VitesseY*NormalY;
  69
  70             minlengthFk=min(minlengthFk,LengthFk/fabs(UN));
  71             minlengthFk=min(minlengthFk,LengthFk/fabs(p.VitesseX));
  72             minlengthFk=min(minlengthFk,LengthFk/fabs(p.VitesseY));
  73
  74             double conc=0.0;
  75             cellCourante=j;
  76             cellAutre=-1;
  77
  78             if (!Fk.isBorder())
  79             {
  80                 int indexC1=Fk.getCellsId()[0];
  81                 int indexC2=Fk.getCellsId()[1];
  82                 /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
  83                 if ( indexC1 == j )
  84                 {
  85                     /* hypothese verifie */
  86                     cellCourante=indexC1;
  87                     cellAutre=indexC2;
  88                 } else if ( indexC2 == j )
  89                 {
  90                     /* hypothese non verifie */
  91                     cellCourante=indexC2;
  92                     cellAutre=indexC1;
  93                 }
  94                 // definir la cellule gauche et droite par le prduit vitesse * normale sortante
  95                 // si u*n>0 : rien a faire sinon inverser la gauche et la droite
  96                 if (UN>1.E-15)
  97                     conc=yField(cellCourante);
  98                 else
  99                     conc=yField(cellAutre);
 100             }else
 101             {
 102                 /* conditions aux limites neumann homogene */
 103                 if (Fk.getGroupName().compare("GAUCHE")==0 || Fk.getGroupName().compare("DROITE")==0)
 104                 {
 105                     if (UN>1.E-15)
 106                         conc=yField(cellCourante);
 107                     else
 108                         conc=0.0;
 109                 }
 110                 /* conditions aux limites periodiques */
 111                 if (Fk.getGroupName().compare("BAS")==0 || Fk.getGroupName().compare("HAUT")==0)
 112                 {
 113                         int indexFP=indexFacesPerio[indexFace];
 114                         /* une autre manière de recuperer l'index de la face periodique */
 115                         //int indexFP=M.getIndexFacePeriodic(indexFace);
 116                         Face Fp=myMesh.getFace(indexFP);
 117                         int indexCp=Fp.getCellsId()[0];
 118                         if (UN>1.E-15)
 119                             conc=yField(cellCourante);
 120                         else
 121                             conc=yField(indexCp);
 122                 }
 123             }
 124             SumF=SumF+UN*LengthFk*conc;
 125           }
 126         dt=p.cfl*minlengthFk/normU;
 127         SumFlux(j)=dt*SumF/Cj.getMeasure();
 128     }
 129 }
 130
 131 void EquationTransport2D(int &iter, double tmin, double & tmax,
 132                          const parameters & p, Field & yField)
 133 {
 134
 135     /*
 136      * MED output of the initial condition at t=0 and iter = 0
 137      */
 138     double time=tmin;
 139     const Mesh& myMesh = yField.getMesh();
 140
 141     yField.setTime(time,iter);
 142     yField.writeMED(p.outFile,true);
 143     yField.writeVTK(p.outFile);
 144     yField.writeCSV(p.outFile);
 145     /* --------------------------------------------- */
 146
 147     /* Time loop */
 148     cout << "Resolution of the transport equation with an UPWIND scheme…" << endl;
 149     int ntmax=3 + iter;
 150     double dt;
 151     IntTab indexFacesPerio=myMesh.getIndexFacePeriodic();
 152     while (iter<ntmax && time <= tmax )
 153     {
 154         Field SumFlux("Fluxes sum",CELLS,myMesh,1) ;
 155         sigma_flux(p,yField,indexFacesPerio,dt,SumFlux);
 156         cout << "-- Iter: " << iter << ", Time: " << time << ", dt: " << dt << endl;
 157
 158         /* Advancing time step */
 159         yField-=SumFlux;
 160
 161         time+=dt;
 162         iter+=1;
 163         // Output every freq iterations
 164         if (iter % p.freqSortie==0)
 165         {
 166             yField.setTime(time,iter);
 167             yField.writeVTK(p.outFile,false);
 168             yField.writeCSV(p.outFile);
 169         }
 170     }
 171
 172     cout << "MED post-treatment of the solution at T=" << time << endl;
 173     yField.setTime(time,iter);
 174     yField.writeMED(p.outFile,false);
 175     yField.writeCSV(p.outFile);
 176
 177     tmax = time;
 178     cout << "End of EquationTransport2D." << endl;
 179 }
 180
 181 void compute_onepass(double tmin, double tmax, const Mesh & M, const parameters & p)
 182 {
 183         int iter = 0;
 184
 185     /* Initial conditions */
 186     cout << "Construction of initial condition…" << endl;
 187     Field f("Y field",CELLS,M,1) ;
 188     conditions_initiales(f);
 189
 190     EquationTransport2D(iter, tmin, tmax, p, f);
 191 }
 192
 193 void compute_twopass(double tmin, double tmax, const Mesh & M, const parameters & p)
 194 {
 195         int iter = 0;
 196         parameters q(p);
 197         double tstep = 0.5 * (tmin + tmax);
 198
 199     /* Initial conditions */
 200     cout << "Construction of initial condition…" << endl;
 201     Field f("Y field",CELLS,M,1) ;
 202     conditions_initiales(f);
 203
 204         q.outFile = p.outFile + "_A";
 205     EquationTransport2D(iter, tmin, tstep, q, f);
 206
 207     Field f2(q.outFile, CELLS, "Y field", iter, 0);
 208
 209         q.outFile = p.outFile + "_B";
 210     EquationTransport2D(iter, tstep, tmax, q, f2);
 211
 212 }
 213
 214 int main()
 215 {
 216     cout << "Resolution of the 2D transport equation:" << endl;
 217     cout << "- DOMAIN: SQUARE" << endl;
 218     cout << "- MESH: TRIANGULAR, GENERATED WITH SALOME" << endl;
 219     cout << "- PERIODIC BC ON TOP AND BOTTOM" << endl;
 220     cout << "- HOMOGENEOUS NEUMANN BC ON LEFT AND RIGHT" << endl;
 221
 222     // donnees du probleme
 223     parameters p;
 224     int iter;
 225
 226     p.cfl=0.4;
 227     p.VitesseX=1.0;
 228     p.VitesseY=1.0;
 229     p.freqSortie=50;
 230     p.outFile="res2D";
 231
 232     cout << "Construction of Cartesian mesh…" << endl;
 233     Mesh M("../../tests/ressources/meshSquare.med");
 234
 235     double tmin = 0;
 236     double tmax = 0.1;
 237
 238     compute_onepass(tmin, tmax, M, p);
 239     compute_twopass(tmin, tmax, M, p);
 240     cout << "CDMATH calculation done." << endl;
 241
 242     return 0;
 243 }