CoreFlows/examples/C/DriftModel_2DInclinedChannelGravityBarriers.cxx

   1 #include "DriftModel.hxx"
   2
   3 using namespace std;
   4
   5 int main(int argc, char** argv)
   6 {
   7         int spaceDim = 2;
   8
   9         //Prepare for the mesh
  10         double xinf=0.0;
  11         double xsup=.6;
  12         double yinf=0.0;
  13         double ysup=2.0;
  14         int nx=3;
  15         int ny=100;
  16         Mesh M(xinf,xsup,nx,yinf,ysup,ny);
  17
  18         //Set the barriers
  19         double xcloison1=xinf+(xsup-xinf)/3;
  20         double xcloison2=xinf+2*(xsup-xinf)/3;
  21         Field barrierField("Barrier Field", FACES, M, 1);
  22         double eps=1e-6;
  23         M.setGroupAtPlan(xsup,0,eps,"wall");
  24         M.setGroupAtPlan(xinf,0,eps,"wall");
  25         M.setGroupAtPlan(ysup,1,eps,"outlet");
  26         M.setGroupAtPlan(yinf,1,eps,"inlet");
  27         double dy=(ysup-yinf)/ny;
  28         int ncloison=3*ny/4;
  29         int i=0;
  30         while( i<= ncloison+1)
  31         {
  32                 M.setGroupAtFaceByCoords(xcloison1,yinf+((ysup-yinf)/4)+(i+0.5)*dy,0,eps,"wall");
  33                 M.setGroupAtFaceByCoords(xcloison2,yinf+((ysup-yinf)/4)+(i+0.5)*dy,0,eps,"wall");
  34                 i++;
  35         }
  36
  37         int nbFaces=M.getNumberOfFaces();
  38         for( i=0;i<nbFaces;i++)
  39         {
  40                 double x=M.getFace(i).x();
  41                 double y=M.getFace(i).y();
  42                 if (((y> yinf+(ysup-yinf)/4) && (abs(x-xcloison1)< eps or abs(x-xcloison2)< eps) ) || abs(x-xinf)< eps || abs(x-xsup)< eps)
  43                         barrierField[i]=1;
  44                 else
  45                         barrierField[i]=0;
  46         }
  47
  48         barrierField.writeVTK("barrierField",true);
  49
  50         // set the limit field for each boundary
  51         double wallVelocityX=0;
  52         double wallVelocityY=0;
  53         double wallTemperature=563;
  54
  55         double inletConcentration=0;
  56         double inletVelocityX=0;
  57         double inletVelocityY=1;
  58         double inletTemperature=563;
  59
  60         double outletPressure=155e5;
  61
  62         // physical constants
  63         vector<double> gravite(spaceDim,0.) ;
  64         gravite[1]=-7;
  65         gravite[0]=7;
  66
  67         DriftModel  myProblem(around155bars600K,spaceDim);
  68         int nVar = myProblem.getNumberOfVariables();
  69
  70         // Prepare for the initial condition
  71         vector<double> VV_Constant(nVar);
  72         // constant vector
  73         VV_Constant[0] = 0;
  74         VV_Constant[1] = 155e5;
  75         VV_Constant[2] = 0;
  76         VV_Constant[3] = 1;
  77         VV_Constant[4] = 563;
  78
  79         //Initial field creation
  80         cout << "Building initial data" << endl;
  81         myProblem.setInitialFieldConstant(M,VV_Constant);
  82
  83         //set the boundary conditions
  84         vector<double>pressure_reference_point(2);
  85         pressure_reference_point[0]=xsup;
  86         pressure_reference_point[1]=ysup;
  87         myProblem.setOutletBoundaryCondition("outlet", outletPressure,pressure_reference_point);
  88         myProblem.setInletBoundaryCondition("inlet", inletTemperature, inletConcentration, inletVelocityX, inletVelocityY);
  89         myProblem.setWallBoundaryCondition("wall", wallTemperature, wallVelocityX, wallVelocityY);
  90
  91         // set physical parameters
  92         myProblem.setGravity(gravite);
  93
  94         // set the numerical method
  95         myProblem.setNumericalScheme(upwind, Explicit);
  96         myProblem.setWellBalancedCorrection(true);
  97         myProblem.setNonLinearFormulation(VFFC);
  98
  99         // name of result file
 100         string fileName = "2DInclinedChannelGravityBarriers";
 101
 102         // computation parameters
 103         unsigned MaxNbOfTimeStep = 3 ;
 104         int freqSave = 1;
 105         double cfl = 0.5;
 106         double maxTime = 500;
 107         double precision = 1e-6;
 108
 109         myProblem.setCFL(cfl);
 110         myProblem.setPrecision(precision);
 111         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
 112         myProblem.setTimeMax(maxTime);
 113         myProblem.setFreqSave(freqSave);
 114         myProblem.setFileName(fileName);
 115         myProblem.usePrimitiveVarsInNewton(true);
 116
 117         // evolution
 118         myProblem.initialize();
 119
 120         bool ok = myProblem.run();
 121         if (ok)
 122                 cout << "Simulation "<<fileName<<" is successful !" << endl;
 123         else
 124                 cout << "Simulation "<<fileName<<"  failed ! " << endl;
 125
 126         cout << "------------ End of calculation !!! -----------" << endl;
 127         myProblem.terminate();
 128
 129         return EXIT_SUCCESS;
 130 }