CoreFlows/examples/C/SinglePhase_1DPorosityJump.cxx

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