CoreFlows/examples/C/DriftModel_1DPorosityJump.cxx

   1 #include "DriftModel.hxx"
   2
   3 using namespace std;
   4
   5 int main(int argc, char** argv)
   6 {
   7         //Preprocessing: mesh and group creation
   8         cout << "Building cartesian mesh" << endl;
   9         double xinf=0.0;
  10         double xsup=4.2;
  11         int nx=100;
  12         Mesh M(xinf,xsup,nx);
  13         double eps=1.E-8;
  14         int spaceDim = M.getSpaceDimension();
  15
  16         //Initial data
  17         double initialConc=0;
  18         double initialVelocityX =1;
  19         double initialTemperature=600;
  20         double initialPressure=155e5;
  21
  22         // physical parameters
  23         Field porosityField("Porosity", CELLS, M, 1);
  24         for(int i=0;i<M.getNumberOfCells();i++){
  25                 double x=M.getCell(i).x();
  26                 if (x> (xsup-xinf)/3 && x< 2*(xsup-xinf)/3)
  27                         porosityField[i]=0.5;
  28                 else
  29                         porosityField[i]=1;
  30         }
  31         porosityField.writeVTK("PorosityField",true);
  32
  33
  34         DriftModel myProblem(around155bars600K,spaceDim);
  35         int nVar = myProblem.getNumberOfVariables();
  36
  37         // Prepare for the initial condition
  38         vector<double> VV_Constant(nVar);
  39         // constant vector
  40         VV_Constant[1] = initialConc;
  41         VV_Constant[1] = initialPressure;
  42         VV_Constant[2] = initialVelocityX;
  43         VV_Constant[3] = initialTemperature;
  44
  45         cout << "Building initial data " << endl;
  46
  47         // generate initial condition
  48         myProblem.setInitialFieldConstant( spaceDim, VV_Constant, xinf, xsup, nx,"Inlet","Outlet");
  49
  50         //set the boundary conditions
  51         myProblem.setInletBoundaryCondition("Inlet",initialTemperature,initialConc,initialVelocityX);
  52         myProblem.setOutletBoundaryCondition("Outlet",initialPressure,vector<double>(1,xsup));
  53
  54         // physical parameters
  55         myProblem.setPorosityField(porosityField);
  56
  57         // set the numerical method
  58         myProblem.setNumericalScheme(upwind, Explicit);
  59         myProblem.setWellBalancedCorrection(true);
  60     myProblem.setNonLinearFormulation(VFFC) ;
  61
  62         // name file save
  63         string fileName = "1DPorosityJumpUpwindWB";
  64
  65
  66         /* set numerical parameters */
  67         unsigned MaxNbOfTimeStep =3;
  68         int freqSave = 1;
  69         double cfl = 0.95;
  70         double maxTime = 5;
  71         double precision = 1e-5;
  72
  73         myProblem.setCFL(cfl);
  74         myProblem.setPrecision(precision);
  75         myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
  76         myProblem.setTimeMax(maxTime);
  77         myProblem.setFreqSave(freqSave);
  78         myProblem.setFileName(fileName);
  79         bool ok;
  80
  81         // evolution
  82         myProblem.initialize();
  83
  84         ok = myProblem.run();
  85         if (ok)
  86                 cout << "Simulation "<<fileName<<" is successful !" << endl;
  87         else
  88                 cout << "Simulation "<<fileName<<"  failed ! " << endl;
  89
  90         cout << "------------ End of calculation -----------" << endl;
  91         myProblem.terminate();
  92
  93         return EXIT_SUCCESS;
  94 }
  95