CoreFlows/examples/C/DriftModel_3DCanalCloison.cxx

   1 #include "DriftModel.hxx"\r
   2 \r
   3 using namespace std;\r
   4 \r
   5 int main(int argc, char** argv)\r
   6 {\r
   7         int spaceDim = 3;\r
   8         // Prepare for the mesh\r
   9         cout << "Building cartesian mesh" << endl;\r
  10         double xinf = 0 ;\r
  11         double xsup=2.0;\r
  12         double yinf=0.0;\r
  13         double ysup=2.0;\r
  14         double zinf=0.0;\r
  15         double zsup=4.0;\r
  16         int nx=10;\r
  17         int ny=nx;\r
  18         int nz=20;\r
  19         double xcloison=(xinf+xsup)/2;\r
  20         double ycloison=(yinf+ysup)/2;\r
  21         double zcloisonmin=1;\r
  22         double zcloisonmax=3;\r
  23         Mesh M(xinf,xsup,nx,yinf,ysup,ny,zinf,zsup,nz);\r
  24         // set the limit field for each boundary\r
  25         double eps=1e-6;\r
  26         M.setGroupAtPlan(xsup,0,eps,"wall");\r
  27         M.setGroupAtPlan(xinf,0,eps,"wall");\r
  28         M.setGroupAtPlan(ysup,1,eps,"wall");\r
  29         M.setGroupAtPlan(yinf,1,eps,"wall");\r
  30         M.setGroupAtPlan(zsup,2,eps,"outlet");\r
  31         M.setGroupAtPlan(zinf,2,eps,"inlet");\r
  32         double dx=(xsup-xinf)/nx;\r
  33         double dy=(ysup-yinf)/ny;\r
  34         double dz=(zsup-zinf)/nz;\r
  35         int ncloison=nz*(zcloisonmax-zcloisonmin)/(zsup-zinf);\r
  36         int i=0 ;\r
  37         int j=0;\r
  38         while( i< ncloison){\r
  39                 while(j< ny){\r
  40                         M.setGroupAtFaceByCoords(xcloison,(j+0.5)*dy,zcloisonmin+(i+0.5)*dz,eps,"wall");\r
  41                         M.setGroupAtFaceByCoords((j+0.5)*dx,ycloison,zcloisonmin+(i+0.5)*dz,eps,"wall");\r
  42                         j=j+1;\r
  43                 }\r
  44                 i=i+1;\r
  45         }\r
  46 \r
  47     // set the limit field for each boundary\r
  48         double wallVelocityX=0;\r
  49         double wallVelocityY=0;\r
  50         double wallVelocityZ=0;\r
  51         double wallTemperature=573;\r
  52         double inletConc=0;\r
  53         double inletVelocityX=0;\r
  54         double inletVelocityY=0;\r
  55         double inletVelocityZ=1;\r
  56         double inletTemperature=563;\r
  57         double outletPressure=155e5;\r
  58 \r
  59     // physical constants\r
  60         vector<double> gravite = vector<double>(spaceDim,0);\r
  61 \r
  62         gravite[0]=0;\r
  63         gravite[1]=0;\r
  64         gravite[2]=-10;\r
  65 \r
  66         double heatPower1=0;\r
  67         double heatPower2=0.25e8;\r
  68         double heatPower3=0.5e8;\r
  69         double heatPower4=1e8;\r
  70 \r
  71         DriftModel myProblem = DriftModel(around155bars600K,spaceDim);\r
  72         int nVar =myProblem.getNumberOfVariables();\r
  73         Field heatPowerField=Field("heatPowerField", CELLS, M, 1);\r
  74 \r
  75         int nbCells=M.getNumberOfCells();\r
  76 \r
  77         for (int i=0;i<nbCells;i++){\r
  78                 double x=M.getCell(i).x();\r
  79                 double y=M.getCell(i).y();\r
  80                 double z=M.getCell(i).z();\r
  81                 if (z> zcloisonmin && z< zcloisonmax)\r
  82                         if (y<ycloison && x<xcloison)\r
  83                                 heatPowerField[i]=heatPower1;\r
  84                         if (y<ycloison && x>xcloison)\r
  85                                 heatPowerField[i]=heatPower2;\r
  86                         if (y>ycloison && x<xcloison)\r
  87                                 heatPowerField[i]=heatPower3;\r
  88                         if (y>ycloison && x>xcloison)\r
  89                                 heatPowerField[i]=heatPower4;\r
  90                 else\r
  91                         heatPowerField[i]=0;\r
  92         }\r
  93         heatPowerField.writeVTK("heatPowerField",true);\r
  94 \r
  95     //Prepare for the initial condition\r
  96         Vector VV_Constant =Vector(nVar);\r
  97 \r
  98         // constant vector\r
  99         VV_Constant[0] = inletConc ;\r
 100         VV_Constant[1] = outletPressure ;\r
 101         VV_Constant[2] = inletVelocityX;\r
 102         VV_Constant[3] = inletVelocityY;\r
 103         VV_Constant[4] = inletVelocityZ;\r
 104         VV_Constant[5] = inletTemperature ;\r
 105 \r
 106     //Initial field creation\r
 107         cout<<"Building initial data " <<endl;\r
 108         myProblem.setInitialFieldConstant(M,VV_Constant);\r
 109 \r
 110     // the boundary conditions\r
 111         myProblem.setOutletBoundaryCondition("outlet", outletPressure);\r
 112         myProblem.setInletBoundaryCondition("inlet", inletTemperature, inletConc, inletVelocityX, inletVelocityY, inletVelocityZ);\r
 113         myProblem.setWallBoundaryCondition("wall", wallTemperature, wallVelocityX, wallVelocityY,wallVelocityZ);\r
 114 \r
 115         // set physical parameters\r
 116         myProblem.setHeatPowerField(heatPowerField);\r
 117         myProblem.setGravity(gravite);\r
 118 \r
 119         // set the numerical method\r
 120         myProblem.setNumericalScheme(upwind, Explicit);\r
 121         myProblem.setWellBalancedCorrection(false);\r
 122 \r
 123         // name file save\r
 124         string fileName = "3DCanalCloison";\r
 125 \r
 126         // parameters calculation\r
 127         unsigned MaxNbOfTimeStep = 3;\r
 128         int freqSave = 1;\r
 129         double cfl = 0.3;\r
 130         double maxTime = 5;\r
 131         double precision = 1e-6;\r
 132 \r
 133         myProblem.setCFL(cfl);\r
 134         myProblem.setPrecision(precision);\r
 135         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);\r
 136         myProblem.setTimeMax(maxTime);\r
 137         myProblem.setFreqSave(freqSave);\r
 138         myProblem.setFileName(fileName);\r
 139         myProblem.setNewtonSolver(precision,20);\r
 140         myProblem.saveVelocity();\r
 141 \r
 142         // evolution\r
 143         myProblem.initialize();\r
 144 \r
 145         bool ok = myProblem.run();\r
 146         if (ok)\r
 147                 cout << "Simulation "<<fileName<<" is successful !" << endl;\r
 148         else\r
 149                 cout << "Simulation "<<fileName<<"  failed ! " << endl;\r
 150 \r
 151         cout << "------------ End of calculation !!! -----------" << endl;\r
 152         myProblem.terminate();\r
 153 \r
 154         return EXIT_SUCCESS;\r
 155 }\r