00001 #include "mainapp.h"
00002 #include "fpclinear.h"
00003 #include <grid/tria.h>
00004 #include <grid/grid_generator.h>
00005 #include "dealbase.h"
00006 #include "boundary_functions_lib.h"
00007 #include "hdf5dealwriter.h"
00008 #include "sequencer.h"
00009 #include "flinear.h"
00010 #include "fquadratic.h"
00011 #include "fchaventmobility.h"
00012 #include "constvelocitymodule.h"
00013 #include "constvectorfunction.h"
00014 #include "fplane.h"
00015 #include "exception.h"
00016 #include "co2gridgenerator.h"
00017 #include "fwellboundarycondition.h"
00018 #include "fchessboard3d.h"
00019 #include "twotracksvectorfunction.h"
00020 #include "facefluxwithgravity.h"
00021 #include "facefluxwithoutgravity.h"
00022 #include "fbucleyleverettgravitymob.h"
00023 #include "fbucleyleveretttotalmob.h"
00024 #include "flinearmobilityproduct.h"
00025 #include "dflinearmobilityproduct.h"
00026 #include "unittests.h"
00027 #include "fconstsphereregion.h"
00028 #include "fdomainunion.h"
00029 #include "fcompoundvector.h"
00030 #include "laxfriedrichsforsystem.h"
00031 #include "transportbase.h"
00032 #include "facefluxwithindependentequations.h"
00033 #include "fluxforco2inj.h"
00034 #include "flashhenrylaw.h"
00035 #include "dummyflash.h"
00036 #include "fdomainunion.h"
00037 #include "upwindforsystem.h"
00038 #include "unittests.h"
00039 #include "russanovforsystem.h"
00040 #include "numericmethods.h"
00041 #include "fwellcondition.h"
00042 #include "fwellcondition.h"
00043 #include <stdarg.h>
00044 #include "godunovmethod.h"
00045 #include "hdf5orthowriter.h"
00046 #include "hdf5orthoreader.h"
00047 #include "stringfunctions.h"
00048 #include "fdomaincomplement.h"
00049 #include "cudaktmethod.h"
00050 #include "ktmethoddbyd.h"
00051 #include "poissonfem.h"
00052 #include "biotfem.h"
00053 #include "flashco2h2o.h"
00054 #include "velreader2d.h"
00055 #include "umfpacksolver.h"
00056 #include "mixedhybridcompositionalsimple.h"
00057 #include "flashco2brine.h"
00058 #include "facefluxco2brine.h"
00059 #include "flashsimpleblackoil.h"
00060 #include "facefluxsimpleblackoil.h"
00061 #include "fconst3d.h"
00062 #include "mixedhybridcompositionalfull.h"
00063 #include "mixedhybridcompositionaltrangenstein1985.h"
00064 #include <typeinfo>
00065 #include "facefluxsimpleblackoilmass.h"
00066 #include "cgsolver.h"
00067 #include "mixedhybridcompositionalsimplepw.h"
00068 #include "fpcsquare.h"
00069 #include "fpcfracture.h"
00070 #include "dfpcfracture.h"
00071 #include "fpcinvfracture.h"
00072 #include <exception>
00073 #include "flashco2brineuncomp.h"
00074 #include "dfpcsquare.h"
00075 #include "solvergpu_agm.h"
00076 #include "mixedhybridincompressible.h"
00077 #include "diffusivestep.h"
00078 #include "compdiffusivestep.h"
00079
00080 #include "blockmatrixmodule.h"
00081 #include "ffractionallinearmobility.h"
00082 #include "fdfractionallinearmobility.h"
00083
00084 #include "finterpolate.h"
00085 #include "fbsinverse.h"
00086 #include "pcinvfunc.h"
00087
00088 #include "fsquaremob.h"
00089 #include "fsquaremobt.h"
00090
00091 #include "stringfunctions.h"
00092 #include "henrylawobiblunt.h"
00093 #include "fwellssource.h"
00094 #include "flashco2brinepw.h"
00095 #include "mumm.h"
00096 #include "biphasicdiff.h"
00097 #include "krpiri1.h"
00098 #include "fmobfromkpermlib.h"
00099
00100 MainApp::MainApp()
00101 {
00102 pOut = &(std::cout);
00103
00104
00105 m_pMesh = NULL;
00106 m_pMB_Mesh = NULL;
00107 m_fTransBC = NULL;
00108 m_fInit = NULL;
00109 m_fK = NULL;
00110 m_pVelFunction = NULL;
00111 m_pDynamicModule = NULL;
00112 m_pTransportModule = NULL;
00113 m_pSequencer = NULL;
00114 m_pFlash = NULL;
00115 m_flux = NULL;
00116 m_fMobT = NULL;
00117 m_fPor = NULL;
00118 m_bPrint=true;
00119 m_blockMatrix = NULL;
00120 }
00121
00122 MainApp::~MainApp()
00123 {
00124
00125
00126
00127
00128
00129
00130 if (m_fTransBC)
00131 delete m_fTransBC;
00132 if (m_fInit)
00133 delete m_fInit;
00134 if (m_fPor)
00135 delete m_fPor;
00136 if (m_fK)
00137 delete m_fK;
00138 if (m_pVelFunction)
00139 delete m_pVelFunction;
00140 if (m_pSequencer)
00141 delete m_pSequencer;
00142 if (m_pFlash)
00143 delete m_pFlash;
00144 if (m_flux)
00145 delete m_flux;
00146 if (m_blockMatrix)
00147 delete m_blockMatrix;
00148 }
00149
00150 Function3D* MainApp::getTransportBC(std::string compName)
00151 {
00152 std::string section="TransportBC";
00153 char keyName[100];
00154 sprintf(keyName,"%s_BOUNDARY_CONDITION",compName.c_str());
00155 print(section,"Boundary Condition for %s: ",compName.c_str());
00156 Function3D *f;
00157 switch(configFile.getInt(keyName))
00158 {
00159 case 0:
00160 {
00161 print(section,"\n\tNo boundary conditions\n");
00162 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00163 configFile.getDouble("DIM_Y"),
00164 configFile.getDouble("DIM_Z"),
00165 INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY);
00166 break;
00167 }
00168 case 1:
00169 {
00170 print(section,"\n\tNo boundary conditions\n");
00171 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00172 configFile.getDouble("DIM_Y"),
00173 configFile.getDouble("DIM_Z"),
00174 INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY);
00175 break;
00176 }
00177 case 2:
00178 {
00179
00180 sprintf(keyName,"%s_BOUNDARY_CONDITION_ARG_1",compName.c_str());
00181 double dA = configFile.getDouble(keyName);
00182 print(section,"\n\t%g at left boundary (%s(0,y,z)=%g)\n",dA,compName.c_str(),dA);
00183 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00184 configFile.getDouble("DIM_Y"),
00185 configFile.getDouble("DIM_Z"),
00186 dA,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY);
00187 break;
00188 }
00189 case 3:
00190 {
00191 print(section,"\n\t1 at right boundary (%s(DIM_X,y,z)=1)\n",compName.c_str());
00192 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00193 configFile.getDouble("DIM_Y"),
00194 configFile.getDouble("DIM_Z"),
00195 INFINITY,1,INFINITY,INFINITY,INFINITY,INFINITY);
00196 break;
00197 }
00198 case 4:
00199 {
00200 print("\n\t1 at front boundary (%s(x,0,z)=1)\n",compName.c_str());
00201 f=new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00202 configFile.getDouble("DIM_Y"),
00203 configFile.getDouble("DIM_Z"),
00204 INFINITY,INFINITY,1,INFINITY,INFINITY,INFINITY);
00205 break;
00206 }
00207 case 5:
00208 {
00209 print(section,"\n\t1 at back boundary (%s(x,DIM_Y,z)=1)\n",compName.c_str());
00210 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00211 configFile.getDouble("DIM_Y"),
00212 configFile.getDouble("DIM_Z"),
00213 INFINITY,INFINITY,INFINITY,1,INFINITY,INFINITY);
00214 break;
00215 }
00216 case 6:
00217 {
00218 print(section,"\n\t1 at bottom boundary (%s(x,y,0)=1)\n",compName.c_str());
00219 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00220 configFile.getDouble("DIM_Y"),
00221 configFile.getDouble("DIM_Z"),
00222 INFINITY,INFINITY,INFINITY,INFINITY,1,INFINITY);
00223 break;
00224
00225 }
00226 case 7:
00227 {
00228 sprintf(keyName,"%s_BOUNDARY_CONDITION_ARG_1",compName.c_str());
00229 double dA = configFile.getDouble(keyName);
00230 print(section,"%g at top boundary (%s(x,y,DIM_Z)=1)\n",dA,compName.c_str());
00231 f = new FCubeBoundaryCondition (configFile.getDouble("DIM_X"),
00232 configFile.getDouble("DIM_Y"),
00233 configFile.getDouble("DIM_Z"),
00234 INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,dA);
00235 break;
00236 }
00237 case 8:
00238 {
00239 char argName[500];
00240 sprintf(argName,"%s_BOUNDARY_CONDITION_ARG_1",compName.c_str());
00241 double height = getMesh().getQ()[2]-getMesh().getP()[2] + getMesh().getGridTolerance();
00242 Point3D C0 = getMesh().getP();
00243 C0[2]-=getMesh().getGridTolerance();
00244 double radius=configFile.getDouble(argName);
00245
00246
00247 print(section,"\n\tFiveSpot boundary condition for %s with radius: %g\n",compName.c_str(),radius);
00248
00249 f = new FCylinderRegion(C0,radius,height,1.0);
00250 break;
00251 }
00252 default:
00253
00254 throw new Exception("\n\tInvalid option for %s_BOUNDARY_CONDITION",compName.c_str());
00255 break;
00256 }
00257 return new FDomainUnion(f, getWellTransportBC(compName),true);
00258 }
00259
00260 Function3D* MainApp::getTransportIC(std::string compName)
00261 {
00262 char keyName[100];
00263 std::string section="IC";
00264 sprintf(keyName,"%s_INITIAL_VALUE",compName.c_str());
00265 print(section,"Initial Condition for %s: ",compName.c_str());
00266 switch (configFile.getInt(keyName))
00267 {
00268 case 1:
00269 print(section,"\n\t%s(x,0) = 0\n",compName.c_str());
00270 return new FPlane(0,0,0,0);
00271 break;
00272 case 2:
00273 {
00274 char arg1[100],arg2[100],arg3[100],arg4[100],arg5[100],arg6[100];
00275 sprintf(arg1,"%s_INITIAL_VALUE_ARG_1",compName.c_str());
00276 sprintf(arg2,"%s_INITIAL_VALUE_ARG_2",compName.c_str());
00277 sprintf(arg3,"%s_INITIAL_VALUE_ARG_3",compName.c_str());
00278 sprintf(arg4,"%s_INITIAL_VALUE_ARG_4",compName.c_str());
00279 sprintf(arg5,"%s_INITIAL_VALUE_ARG_5",compName.c_str());
00280 sprintf(arg6,"%s_INITIAL_VALUE_ARG_6",compName.c_str());
00281
00282 print(section,"\n\t%s equal to %g in a sphere region of radius %g centered in <%g, %g, %g>\n\t and %g everywhere else.",
00283 keyName,
00284 configFile.getDouble(arg1),
00285 configFile.getDouble(arg3),
00286 configFile.getDouble(arg4),
00287 configFile.getDouble(arg5),
00288 configFile.getDouble(arg6),
00289 configFile.getDouble(arg2));
00290
00291
00292 return new FConstSphereRegion(configFile.getDouble(arg1),
00293 configFile.getDouble(arg2),
00294 configFile.getDouble(arg3),
00295 configFile.getDouble(arg4),
00296 configFile.getDouble(arg5),
00297 configFile.getDouble(arg6));
00298 break;
00299 }
00300 case 3:
00301 {
00302 char keyFileName[1000];
00303 sprintf(keyFileName,"%s_INITIAL_VALUE_ARG_1",compName.c_str());
00304 Point3D p0(0,0,0);
00305 Point3D p1(configFile.getDouble("DIM_X"),
00306 configFile.getDouble("DIM_Y"),
00307 configFile.getDouble("DIM_Z"));
00308 FChessBoard3D *f = new FChessBoard3D(p0,p1,getMesh().getP(),getMesh().getQ(),configFile.getString(keyFileName));
00309 print(section,"Chess Board function.\n\tfile: \"%s\"\n\tDimensions: %d %d %d\n",
00310 configFile.getString(keyFileName).c_str(),f->getDimX(),f->getDimY(),f->getDimZ());
00311 return f;
00312 }
00313 case 4:
00314 {
00315 char Arg1[1000];
00316 sprintf(Arg1,"%s_INITIAL_VALUE_ARG_1",compName.c_str());
00317 double dd=configFile.getDouble(Arg1);
00318 print(section,"\n\t%s(x,0) = %g\n",compName.c_str(),dd);
00319 return new FPlane(0.0,0.0,0.0,dd);
00320 }
00321
00322 default:
00323
00324 throw new Exception("Invalid option for %s_INITIAL_VALUE",compName.c_str());
00325 return NULL;
00326 break;
00327 }
00328 return NULL;
00329 }
00330
00334 Function3D* MainApp::getPressureIC()
00335 {
00336
00337 std::string section="PressureIC";
00338
00339 print(section,"Initial Condition for Pressure: ");
00340
00341 switch (configFile.getInt("PRESSURE_INITIAL_VALUE"))
00342 {
00343 case 1:
00344 {
00345 double p = configFile.getDouble("PRESSURE_INITIAL_ARG_1");
00346 print(section,"\n\tp(x,0) = %g N/m^2\n",p);
00347 m_fInitPressure= new FPlane(0,0,0,p);
00348 return m_fInitPressure;
00349 break;
00350 }
00351 default:
00352 throw new Exception("Invalid option for PRESSURE_INITIAL_VALUE");
00353 return NULL;
00354 break;
00355 }
00356 return NULL;
00357 }
00358
00359 Sequencer* MainApp::getSequencer()
00360 {
00361
00362 return new Sequencer(getMonitorWells());
00363 }
00364
00365 void MainApp::setTriangulation(Triangulation<3>& tria)
00366 {
00367 Point<3> P1(0,0,0);
00368 Point<3> P2(configFile.getDouble("DIM_X"),configFile.getDouble("DIM_Y"),configFile.getDouble("DIM_Z"));
00369
00370 std::vector<unsigned> v(3);
00371 v[0]=(unsigned) configFile.getDouble("NUM_ELEMS_X");
00372 v[1]=(unsigned) configFile.getDouble("NUM_ELEMS_Y");
00373 v[2]=(unsigned) configFile.getDouble("NUM_ELEMS_Z");
00374
00375 VecWellInfo wells = getWells();
00376
00377
00378 CO2GridGenerator::subdivided_hyper_rectangle_with_hole(tria,v,P1,P2,wells);
00379 }
00380
00381 OrthoMesh& MainApp::getMesh()
00382 {
00383 if (!m_pMesh)
00384 {
00385 Point<3> P1(0,0,0);
00386 Point<3> P2(configFile.getDouble("DIM_X"),configFile.getDouble("DIM_Y"),configFile.getDouble("DIM_Z"));
00387
00388 VecWellInfo wells = getWells();
00389
00390
00391 VecWellInfo::iterator it = wells.begin();
00392 while (it!=wells.end())
00393 {
00394 if (it->getBCType() == WellInfo::PRESSURE || it->getBCType() == WellInfo::FLUX)
00395 it++;
00396 else
00397 it=wells.erase(it);
00398 }
00399
00400 m_pMesh = new OrthoMesh(P1,P2,
00401 configFile.getInt("NUM_ELEMS_X"),
00402 configFile.getInt("NUM_ELEMS_Y"),
00403 configFile.getInt("NUM_ELEMS_Z"),wells);
00404 printTriangulation();
00405 }
00406 return *m_pMesh;
00407 }
00408
00409 OrthoMesh& MainApp::getMB_Mesh()
00410 {
00411 if (!m_pMB_Mesh)
00412 {
00413 Point<3> P1(0,0,0);
00414 Point<3> P2(configFile.getDouble("MB_DIM_X"),configFile.getDouble("MB_DIM_Y"),configFile.getDouble("MB_DIM_Z"));
00415
00416 m_pMB_Mesh = new OrthoMesh(P1,P2,
00417 configFile.getInt("MB_NUM_ELEMS_X"),
00418 configFile.getInt("MB_NUM_ELEMS_Y"),
00419 configFile.getInt("MB_NUM_ELEMS_Z"));
00420 HDF5OrthoWriter::getHDF5OrthoWriter().writeMesh("MBmesh",*m_pMB_Mesh);
00421
00422 }
00423 return *m_pMB_Mesh;
00424 }
00425
00426
00427 TransportBase* MainApp::getTransportModule()
00428 {
00429 std::string section="Transport";
00430 if (m_pTransportModule)
00431 return m_pTransportModule;
00432
00433
00434
00435
00436 switch (configFile.getInt("TRANSPORT_MODULE"))
00437 {
00438 case 1:
00439 {
00440 print(section,"Method: Lax Friedrichs\nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00441 OrthoMesh &mesh = getMesh();
00442 m_flux=getFluxFunction();
00443 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00444 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00445 m_pTransportModule= new LaxFriedrichsForSystem(mesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,getTransportFixedConditions(),configFile.getDouble("MAX_CFL"));
00446 break;
00447 }
00448 case 2:
00449 {
00450 print(section,"Method: Upwind \nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00451 OrthoMesh &mesh = getMesh();
00452 m_flux=getFluxFunction();
00453 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00454 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00455 DiffusiveStep *diffStp = getDiffusiveStep();
00456 m_pTransportModule= new UpwindForSystem(mesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,getTransportFixedConditions(),configFile.getDouble("MAX_CFL"),diffStp);
00457 break;
00458 }
00459 case 3:
00460 {
00461 print(section,"Method: Russanov \nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00462 OrthoMesh &mesh = getMesh();
00463 m_flux=getFluxFunction();
00464 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00465 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00466 m_pTransportModule = new RussanovForSystem(mesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,getTransportFixedConditions(),configFile.getDouble("MAX_CFL"));
00467 break;
00468 }
00469 case 4:
00470 print(section,"Method: Godunov\nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00471 m_flux=getFluxFunction();
00472 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00473 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00474 m_pTransportModule= new GodunovMethod(*m_pMesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,getTransportFixedConditions(),configFile.getDouble("MAX_CFL"));
00475 break;
00476 case 5:
00477 print(section,"Method: Kurganov-Tadmor (KT) 2nd Order\nDimension by Dymension\nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00478 m_flux=getFluxFunction();
00479 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00480 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00481 m_pTransportModule= new KTMethodDByD(*m_pMesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,*m_flux,getTransportFixedConditions(),configFile.getDouble("MAX_CFL"));
00482 break;
00483 case 6:
00484 {
00485 print(section,"Method: Kurganov-Tadmor (KT) 2nd Order (by The GPU Master Rahunassalan)\nCFL: %f\n",configFile.getDouble("MAX_CFL"));
00486 m_flux=getFluxFunction();
00487 m_fTransBC = getTransportForSystemBC(m_flux->numComponents());
00488 m_fInit = getTransportForSystemIC(m_flux->numComponents());
00489 double v1 = configFile.getDouble("S1_VISCOSITY");
00490 double v2 = configFile.getDouble("S2_VISCOSITY");
00491 double theta = configFile.getDouble("KT_THETA");
00492 double rs1 = configFile.getDouble("RESIDUAL_S1");
00493 double rs2 = configFile.getDouble("RESIDUAL_S2");
00494 print(section,"Mobility Funcion: Bucley Leverret with viscosities %g,%g\n",v1,v2);
00495 m_pTransportModule= new CUDAKTMethod(*m_pMesh,*m_fInit,getPorosityAtCells(),*m_fTransBC,v1,v2,theta,rs1,rs2,configFile.getDouble("MAX_CFL"));
00496 break;
00497 }
00498 default:
00499
00500 throw new Exception("Invalid option for TRANSPORT_MODULE selection");
00501 return NULL;
00502
00503 }
00504 return m_pTransportModule;
00505 }
00506
00507
00508 FlashBase* MainApp::getFlash()
00509 {
00510 std::string section="Flash";
00511
00512 if (m_pFlash!=NULL)
00513 return m_pFlash;
00514 switch(configFile.getInt("FLASH_MODULE"))
00515 {
00516 case 1:
00517 {
00518 print(section,"Dummy Flash that does not do any computation\n");
00519 unsigned n_phases = configFile.getInt("DUMMY_FLASH_NUM_PHASES");
00520 unsigned n_comps = configFile.getInt("DUMMY_FLASH_NUM_COMPONENTS");
00521 VecDouble viscs(n_phases);
00522 VecDouble mass(n_comps);
00523 char name[2000];
00524 print(section,"Num Components %d\n",n_comps);
00525 print(section,"Num Phases %d\n",n_phases);
00526 print(section,"Viscosities: ");
00527 for (unsigned i=1;i<=n_phases;i++)
00528 {
00529 sprintf(name,"S%d_VISCOSITY",i);
00530 viscs(i-1)=configFile.getDouble(name);
00531 print(section,"\n\tPhase %d = %g ",i,viscs(i-1));
00532 }
00533 print(section,"\nDensities: ");
00534 for (unsigned i=1;i<=n_comps;i++)
00535 {
00536 sprintf(name,"S%d_DENSITY",i);
00537 mass(i-1)=configFile.getDouble(name);
00538 print(section,"\n\tPhase %d = %g ",i,mass(i-1));
00539 }
00540
00541
00542 print(section,"\n");
00543 m_pFlash= new DummyFlash(n_phases,n_comps,viscs,mass);
00544 break;
00545 }
00546 case 2:
00547 {
00548 print(section,"Flash for the Henry Law computation of dissolution of gas into liquid\n");
00549 print(section,"\tXd = P/Kh exp(-vP/(R*T))\n");
00550 print(section,"\tKd = Xd mL/(1-Xd(1-mL/mc))\n");
00551 print(section,"\tKd is the maximum concentration of gas moles in the liquid\n");
00552 print(section,"\tHenry Law states that the dissolution is instantly and the\n\t\
00553 concentration of the dissolved gas is always Kd except for the \n\t\
00554 case there is no enough gas to dissolve. In this case , all the gas is\n\t\
00555 dissolved in the liquid\n");
00556 print(section,"\tHENRY_COEFFICIENT Kh = %g\n",configFile.getDouble("HENRY_COEFFICIENT"));
00557 print(section,"\tAPPARENT CO2 MOLAR VOLUME v = %g\n",configFile.getDouble("APPARENT_CO2_MOLAR_VOLUME"));
00558 print(section,"\tGAS UNIVERSAL CONSTANT R = = %g\n",configFile.getDouble("R"));
00559 print(section,"\tCONSTANT TEMPERATURE T = %g\n",configFile.getDouble("TEMPERATURE"));
00560 print(section,"\tGAS MOLAR VOLUME mc = %g\n",configFile.getDouble("GAS_MOLAR_VOLUME"));
00561 print(section,"\tLIQUID MOLAR VOLUME mL = %g\n",configFile.getDouble("LIQUID_MOLAR_VOLUME"));
00562 m_pFlash=new FlashHenryLaw(getMesh(),
00563 configFile.getDouble("HENRY_COEFFICIENT"),
00564 configFile.getDouble("APPARENT_CO2_MOLAR_VOLUME"),
00565 configFile.getDouble("R"),
00566 configFile.getDouble("TEMPERATURE"),
00567 configFile.getDouble("GAS_MOLAR_VOLUME"),
00568 configFile.getDouble("LIQUID_MOLAR_VOLUME"),
00569 configFile.getDouble("S1_VISCOSITY"),
00570 configFile.getDouble("S2_VISCOSITY"));
00571
00572 break;
00573 }
00574 case 3:
00575 {
00576 print(section,"Flash for the Henry Law computation of dissolution of gas into liquid\n");
00577 print(section,"designed for the Obi Blunt model that is formulated in terms of Saturation\n\
00578 of CO2 and concentration of CO2 in water");
00579 print(section,"\tXd = P/Kh exp(-vP/(R*T))\n");
00580 print(section,"\tKd = Xd mL/(1-Xd(1-mL/mc))\n");
00581 print(section,"\tKd is the maximum concentration of gas moles in the liquid\n");
00582 print(section,"\tHenry Law states that the dissolution is instantly and the\n\t\
00583 concentration of the dissolved gas is always Kd except for the \n\t\
00584 case there is no enough gas to dissolve. In this case , all the gas is\n\t\
00585 dissolved in the liquid\n");
00586 print(section,"\tHENRY_COEFFICIENT Kh = %g\n",configFile.getDouble("HENRY_COEFFICIENT"));
00587 print(section,"\tAPPARENT CO2 MOLAR VOLUME v = %g\n",configFile.getDouble("APPARENT_CO2_MOLAR_VOLUME"));
00588 print(section,"\tGAS UNIVERSAL CONSTANT R = = %g\n",configFile.getDouble("R"));
00589 print(section,"\tCONSTANT TEMPERATURE T = %g\n",configFile.getDouble("TEMPERATURE"));
00590 print(section,"\tGAS MOLAR VOLUME mc = %g\n",configFile.getDouble("GAS_MOLAR_VOLUME"));
00591 print(section,"\tLIQUID MOLAR VOLUME mL = %g\n",configFile.getDouble("LIQUID_MOLAR_VOLUME"));
00592 m_pFlash= new HenryLawObiBlunt(getMesh(),
00593 configFile.getDouble("HENRY_COEFFICIENT"),
00594 configFile.getDouble("APPARENT_CO2_MOLAR_VOLUME"),
00595 configFile.getDouble("R"),
00596 configFile.getDouble("TEMPERATURE"),
00597 configFile.getDouble("GAS_MOLAR_VOLUME"),
00598 configFile.getDouble("LIQUID_MOLAR_VOLUME"));
00599 break;
00600
00601
00602 }
00603 case 4:
00604 {
00605
00606 double T = configFile.getDouble("TEMPERATURE");
00607 print(section,"\
00608 Flash Equilibrium Method for CO2 Storage in Brine Aquiffers. \n\
00609 Authors: Allan Leal and Mohammad Peri, 2010\n\
00610 TEMPERATURE: %g Kelvin\n\
00611 ",T);
00612
00613
00614
00615 m_pFlash= new FlashCO2Brine(getMesh(),T);
00616 break;
00617 }
00618 case 5:
00619 {
00620 print(section,"\
00621 Simple Black Oil Model for OIL and GAS where only gas\n \
00622 can dissolve in oil. The units for this flash is\n\
00623 distance ft\n\
00624 time days\n\
00625 Pressure PSI\n\
00626 \n\
00627 This flash override the viscosities values given to\n\
00628 the mobility functions such that they depend on\n\
00629 the pressure. So the viscosity values given to the\n\
00630 mobility functions are not relevant.\n\
00631 \"John B Bell, Gregory R. Shubin, John Trangestein,\n\
00632 METHOD FOR REDUCING NUMERICAL DISPERSION IN TWO\n\
00633 PHASE BLACK OIL RESERVOIR SIMULATION, Journal of\n\
00634 Computational Physics 65, 71-106 1986\"\n");
00635 m_pFlash= new FlashSimpleBlackOil(getMesh());
00636 }
00637 case 6:
00638 {
00639
00640 double T = configFile.getDouble("TEMPERATURE");
00641 print(section,"\
00642 Uncompressible Modification of the \n\
00643 Flash Equilibrium Method for CO2 Storage in Brine Aquiffers. \n\
00644 Authors: Allan Leal and Mohammad Peri, 2010\n\
00645 TEMPERATURE: %g Kelvin\n\
00646 Where all the components have the same molar volume = 1\n\
00647 nomatter wich phase. In this way we turn of the compressibility\n\
00648 without altering the dissolution.\n\
00649 ",T);
00650 m_pFlash= new FlashCO2BrineUncomp(getMesh(),T);
00651 }
00652
00653 case 7:
00654 {
00655 double T = configFile.getDouble("TEMPERATURE");
00656 print(section,"\
00657 Flash Equilibrium Method for CO2 Storage in Brine Aquiffers. \n\
00658 Authors: Allan Leal and Mohammad Peri, 2010\n\
00659 TEMPERATURE: %g Kelvin\n\
00660 This is an extension of this flash in order to incorporate capillary\n\
00661 pressure. The pressure passed on this module is actually the water pressure,\n\
00662 that is converted internally to be Sw*pw + So*po"
00663 ,T);
00664
00665
00666
00667 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
00668 getPcFunction(&fPc,&dfPc,&fPcInv);
00669
00670 m_pFlash= new FlashCO2BrinePw(getMesh(),T,dynamic_cast<Function1D&>(*fPc));
00671 break;
00672
00673
00674 }
00675 default:
00676 throw new Exception("Invalid option for key FLASH_MODULE");
00677 }
00678
00679 return m_pFlash;
00680 }
00681
00682 void MainApp::execute(std::string fileName,bool bunittests)
00683 {
00684
00685 try {
00686
00687 configFile.readFile(fileName);
00688
00689
00690 if (bunittests)
00691 {
00692 UnitTests tests(*this);
00693 tests.execute();
00694 printBeginSection("END UNITS TESTS");
00695 return;
00696 }
00697
00698
00699 HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
00700 hdf5.setOutputFile(getHDF5FileName());
00701 print("HDF5","HDF5 Ouput: %s",getHDF5FileName().c_str());
00702 hdf5.writeMesh("tria0",getMesh());
00703 hdf5.setVariable("Time",0);
00704
00705 m_pSequencer = getSequencer();
00706
00707 switch (configFile.getInt("SEQUENCE_ALGORITHM"))
00708 {
00709 case 0:
00710 print("Sequencer","Dynamic Module Test\n");
00711 printLog(std::cout);
00712 m_pTransportModule = getTransportModule();
00713 m_pDynamicModule = getDynamicModule();
00714 m_pSequencer->testDynamicModule(*m_pDynamicModule,*m_pTransportModule);
00715 break;
00716 case 1:
00717 {
00718 m_pTransportModule = getTransportModule();
00719 m_pDynamicModule = getDynamicModule();
00720
00721 print("Sequencer","Method: Transport Test\nEnd Time: %g\nOutputs: %g\n\n\n",
00722 configFile.getDouble("END_TIME"),
00723 configFile.getDouble("N_OUTPUTS"));
00724 printLog(std::cout);
00725
00726 m_pSequencer->testTransport(*m_pDynamicModule,*m_pTransportModule,
00727 configFile.getDouble("END_TIME"),
00728 configFile.getDouble("N_OUTPUTS"));
00729 break;
00730 }
00731 case 2:
00732 {
00733 m_pTransportModule = getTransportModule();
00734 m_pDynamicModule = getDynamicModule();
00735 m_pDiff = getDiffusiveStep();
00736 if (m_pDiff)
00737 {
00738 m_pDiff->setTransport(*m_pTransportModule);
00739 m_pDiff->setDynamic(*m_pDynamicModule);
00740 }
00741 int nDyn = configFile.getInt("N_DYN_ITERATIONS");
00742 int nOuts = configFile.getInt("N_OUTPUTS");
00743 int nOutsPerDyn = nOuts/nDyn;
00744 if (nOutsPerDyn == 0)
00745 nOutsPerDyn = 1;
00746 print("Sequencer","Method: Staggered Iteration\nEnd Time: %g\nDynamic Modules Iterations: %d\nOutputs Per Dyn Iteration: %d\n",
00747 configFile.getDouble("END_TIME"),
00748 configFile.getInt("N_DYN_ITERATIONS"),
00749 nOutsPerDyn);
00750 printLog(std::cout);
00751 m_pSequencer->alternateIteration(*m_pDynamicModule,*m_pTransportModule,m_pDiff,
00752 configFile.getDouble("END_TIME"),
00753 configFile.getInt("N_DYN_ITERATIONS"),
00754 nOuts);
00755 break;
00756 }
00757 case 3:
00758 {
00759 m_pFlash = getFlash();
00760 m_pTransportModule = getTransportModule();
00761 m_pDynamicModule = getDynamicModule();
00762
00763 m_pTransportModule->setFlash(dynamic_cast<FlashCompositional*>(getFlash()));
00764 m_pFlash->setTransport(*m_pTransportModule);
00765 m_pFlash->setDynamic(*m_pDynamicModule);
00766
00767 m_pDiff = getDiffusiveStep();
00768 if (m_pDiff != NULL)
00769 {
00770 m_pDiff->setTransport(*m_pTransportModule);
00771 m_pDiff->setDynamic(*m_pDynamicModule);
00772 }
00773
00774
00775
00776
00777
00778
00779 if (dynamic_cast<FlashCompositional*>(getFlash()))
00780 {
00781 print("Sequencer",
00782 "Adjusting Initial Conditions for Compressible Model\n");
00783 adjustInitialConditionForCompressibleModel(*m_pDynamicModule,dynamic_cast<ConservativeMethodForSystem&>(*m_pTransportModule),*getFlashCompositional());
00784 }
00785
00786
00787
00788
00789 print("Sequencer","Method: Staggered Iteration with Flash calculation\n");
00790 print("Sequencer","End Time: %g\n" ,configFile.getDouble("END_TIME"));
00791 print("Sequencer","Dynamic Iterations: %g\n",configFile.getDouble("N_DYN_ITERATIONS"));
00792 ConservativeMethodForSystem& transportSystem = dynamic_cast<ConservativeMethodForSystem&>(*m_pTransportModule);
00793
00794
00795 if (configFile.isDefined("TRANSPORT_PER_DYNAMIC_STEPS"))
00796 {
00797 CompressibleDynamic *pCompDynamic;
00798 pCompDynamic = dynamic_cast<CompressibleDynamic*>(&(*m_pDynamicModule));
00799 if (!pCompDynamic)
00800 {
00801 throw new Exception("The proportion control sequence alghorithm demands a module that implements CompressibleDynamicBase\n");
00802 }
00803 print("Sequencer","Max Transport Steps per Cycle: %g\n",configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS"));
00804 printLog(std::cout);
00805 m_pSequencer->alternateIterationProportionControl(*pCompDynamic,
00806 transportSystem,
00807 *m_pFlash,m_pDiff,
00808 configFile.getDouble("END_TIME"),
00809 configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS"),
00810 configFile.getDouble("TRANSPORT_PER_DYNAMIC_STEPS_TOLERANCE"),
00811 getDynamicTimeStep(),
00812 configFile.getDouble("N_OUTPUTS"));
00813
00814 }
00815 else
00816 {
00817
00818 printLog(std::cout);
00819 m_pSequencer->alternateIteration(*m_pDynamicModule,
00820 transportSystem,
00821 *m_pFlash,
00822 m_pDiff,
00823 configFile.getDouble("END_TIME"),
00824 configFile.getDouble("N_DYN_ITERATIONS"),
00825 configFile.getDouble("N_OUTPUTS"));
00826 }
00827 }
00828 break;
00829 case 4:
00830 {
00831 m_pTransportModule = getTransportModule();
00832 m_pDiff = getDiffusiveStep();
00833 m_pDiff->setTransport(*m_pTransportModule);
00834
00835 int nDyn = configFile.getInt("N_DYN_ITERATIONS");
00836 unsigned nOuts = configFile.getInt("N_OUTPUTS");
00837 int nOutsPerDyn = nOuts/nDyn;
00838 if (nOutsPerDyn == 0)
00839 nOutsPerDyn = 1;
00840 printLog(std::cout);
00841 m_pSequencer->diffusiveStepTest(*m_pDiff, configFile.getInt("N_DYN_ITERATIONS"),
00842 configFile.getDouble("END_TIME"),
00843 nOuts);
00844
00845 break;
00846
00847
00848 }
00849 default:
00850 throw new Exception("Wrong option for the SEQUENCE_ALGORITHM key\n");
00851 return;
00852 }
00853 printf("\n\nSimulation Done\n");
00854 hdf5.close();
00855 return;
00856
00857 }
00858 catch(Exception *e)
00859 {
00860 printf("Error");
00861 printBeginSection("ERROR MESSAGE");
00862 printf("%s\n\nQuitting.....\n",e->getError().c_str());
00863 abort();
00864 }
00865
00866 }
00867
00872 DynamicBase* MainApp::getDynamicModule()
00873 {
00874 std::string section="Dynamic";
00875 std::string filename;
00876 if (m_pDynamicModule)
00877 return m_pDynamicModule;
00878
00879 int dynMod = configFile.getInt("DYNAMIC_MODULE");
00880 switch (dynMod)
00881 {
00882 case 1:
00883 print(section,"Module that gives an constant velocity field.\nIt is used mainly for test purposes\n");
00884 m_pDynamicModule= this->getConstVelocityDynamicModule();
00885 break;
00886 case 2:
00887 {
00888 print(section,"Modules that read the 2D velocity fluid from a file and project it in 3D ");
00889 int eleZ = configFile.getInt("NUM_ELEMS_Z");
00890 if(eleZ != 1)
00891 {
00892 throw new Exception("The number of element in Z direction must be 1");
00893 exit(1);
00894 }
00895 filename= configFile.getString("VELOCITY_FILE");
00896 print("File name: %s\n",filename.c_str());
00897 m_pDynamicModule= new VelReader2D(filename,getMesh());
00898 break;
00899
00900 }
00901 case 3:
00902 {
00903 print(section,"Module that uses the mixed hybrid finite element method\nfor compositional models\n");
00904 getPBoundaryCondition(m_fDP,m_fFluxBC);
00905 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
00906 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
00907 VecDouble &vPor=getPorosityAtCells();
00908 m_fInitPressure=getPressureIC();
00909 double grav=configFile.getDouble("GRAVITY");
00910 double dt =getDynamicTimeStep();
00911 m_pSolver =getLinearSolver();
00912 print(section,"GRAVITY: %g\n",grav);
00913 print(section,"TIME STEP: %g\n",dt);
00914
00915 m_pDynamicModule= new MixedHybridCompositionalSimple(getMesh(),*m_fInitPressure,*m_fFluxBC,*m_fDP,getPermeabilityAtCells(),vPor, dynamic_cast<Function1D&>(*fMobT),dynamic_cast<Function1D&>(*fFracF),grav,dt,*getFlashCompositional(),*m_pSolver,configFile.getInt("DYNAMIC_DEBUG_LEVEL"));
00916 break;
00917 }
00918 case 4:
00919 {
00920
00921 print(section,"Biot Problem Test");
00922 getPBoundaryCondition(m_fDP,m_fFluxBC);
00923 VecDouble &vK = getPermeabilityAtCells();
00924 m_pSolver =getLinearSolver();
00925 m_pDynamicModule= new PoissonFEM(getMesh(),vK,*m_fDP,*m_fFluxBC,*getLinearSolver() );
00926 break;
00927 }
00928 case 5:
00929 {
00930 print(section,"Module that uses the mixed hybrid finite element method\nfor full compositional models\n");
00931 getPBoundaryCondition(m_fDP,m_fFluxBC);
00932 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
00933 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
00934 VecDouble &vPor=getPorosityAtCells();
00935 m_fInitPressure=getPressureIC();
00936 double grav=configFile.getDouble("GRAVITY");
00937 double dt =getDynamicTimeStep();
00938 m_pSolver =getLinearSolver();
00939 print(section,"GRAVITY: %g\n",grav);
00940 print(section,"TIME STEP: %g\n",dt);
00941
00942 m_pDynamicModule= new MixedHybridCompositionalFull(getMesh(),*m_fInitPressure,*m_fFluxBC,*m_fDP,getPermeabilityAtCells(),vPor, dynamic_cast<Function1D&>(*fMobT),dynamic_cast<Function1D&>(*fFracF),grav,dt,*getFlashCompositional(),*m_pSolver,configFile.getInt("DYNAMIC_DEBUG_LEVEL"));
00943 break;
00944 }
00945 case 6:
00946 {
00947 print(section,"Module that solve the hybrid finite element for biphase incompressible flow\n");
00948 double grav = configFile.getDouble("GRAVITY");
00949 VecDouble densities(2);
00950 densities(0)=configFile.getDouble("S1_DENSITY");
00951 densities(1)=configFile.getDouble("S2_DENSITY");
00952
00953 print(section,"Gravity: %g\n",grav);
00954 print(section,"Phases Densities:\n\tPho1 = %g\n\tPho2 = %g\n\t",densities(0),densities(1));
00955
00956
00957 getPBoundaryCondition(m_fDP,m_fFluxBC);
00958 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
00959 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
00960 m_pSolver = getLinearSolver();
00961
00962
00963
00964
00965 m_pDynamicModule= new MixedHybridIncompressible(getMesh(),*m_fFluxBC,*m_fDP,
00966 getPermeabilityAtCells(),
00967 *fMobT,*fFracF,
00968 grav,densities,
00969 *m_pSolver,
00970 configFile.getInt("DYNAMIC_DEBUG_LEVEL"));
00971 break;
00972 }
00973 case 7:
00974 {
00975 print(section,"Module that uses the mixed hybrid finite element method considering cappillar pressure and using Pw as main variable\n");
00976 getPBoundaryCondition(m_fDP,m_fFluxBC);
00977
00978 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
00979
00980 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
00981 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
00982 getPcFunction(&fPc,&dfPc,&fPcInv);
00983 VecDouble &vPor=getPorosityAtCells();
00984 m_fInitPressure=getPressureIC();
00985 enablePrint(false);
00986 Function3D *fTBC = getTransportForSystemBC(getFlash()->numComponents());
00987 enablePrint(true);
00988 double grav=configFile.getDouble("GRAVITY");
00989 double dt =getDynamicTimeStep();
00990 m_pSolver =getLinearSolver();
00991 print(section,"GRAVITY: %g\n",grav);
00992 print(section,"TIME STEP: %g\n",dt);
00993
00994
00995 m_pDynamicModule= new MixedHybridCompositionalSimplePw(getMesh(),*m_fInitPressure,*m_fFluxBC,*m_fDP,*fTBC,getPermeabilityAtCells(),vPor, *fMobT,*fFracF,*fPc,grav,dt,*getFlashCompositional(),*m_pSolver,configFile.getInt("DYNAMIC_DEBUG_LEVEL"));
00996 break;
00997
00998 }
00999 case 8:
01000 {
01001 m_pDynamicModule= &getBlockMatrix();
01002 break;
01003 }
01004 case 9:
01005 {
01006 print(section,"This is the MUMM method for 2D\nBe sure the mesh has just one element in Z\n");
01007
01008 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
01009 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
01010
01011 m_pDynamicModule= new MUMM(getMesh(),getPermeabilityAtCells(),dynamic_cast<Function1D&>(*fMobT));
01012 break;
01013 }
01014
01015 case 11:
01016 {
01017 print(section,"Elliptic equation solved in terms of the pressure and Petrov-Galerkin\npost processing\
01018 of the velocity fields\n");
01019 getPBoundaryCondition(m_fDP,m_fFluxBC);
01020
01021 m_pDynamicModule= new PoissonFEM(getMesh(),getPermeabilityAtCells(),
01022 *m_fDP,*m_fFluxBC,
01023 *getLinearSolver(),
01024 configFile.getUnsigned("DYNAMIC_DEBUG_LEVEL"));
01025
01026
01027
01028
01029 break;
01030
01031 }
01032 case 12:
01033 {
01034 print(section,"Biot Problem solved with UMFPACK\n");
01035 getPBoundaryCondition(m_fDP,m_fFluxBC);
01036 getUBoundaryCondition(m_fDU,m_fTensionBC);
01037 double dYoung = getYoung();
01038 double dPoisson = getPoisson();
01039 double dt = getDynamicTimeStep();
01040
01041
01042 m_pDynamicModule= new BiotFEM(getMesh(),dt,getPermeabilityAtCells(),dYoung,dPoisson,*m_fDU,*m_fTensionBC);
01043 break;
01044
01045 }
01046 case 13:
01047 {
01048 print(section,"Module that uses the mixed hybrid finite element method for compositional models\n");
01049 print(section,"explained in Bell Trangenstein in 1985\n");
01050 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
01051 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
01052 getPBoundaryCondition(m_fDP,m_fFluxBC);
01053 VecDouble &vPor=getPorosityAtCells();
01054 m_fInitPressure=getPressureIC();
01055 double grav=configFile.getDouble("GRAVITY");
01056 double dt =getDynamicTimeStep();
01057 m_pSolver =getLinearSolver();
01058 print(section,"GRAVITY: %g\n",grav);
01059 print(section,"TIME STEP: %g\n",dt);
01060
01061 m_pDynamicModule= new MixedHybridCompositionalTrangenstein1985(getMesh(),*m_fInitPressure,*m_fFluxBC,*m_fDP,getPermeabilityAtCells(),vPor, dynamic_cast<Function1D&>(*fMobT),dynamic_cast<Function1D&>(*fFracF),grav,dt,*getFlashCompositional(),*m_pSolver,configFile.getInt("DYNAMIC_DEBUG_LEVEL"));
01062 break;
01063 }
01064
01065 default:
01066 {
01067 if (dynMod >= 2 && dynMod <= 5)
01068 {
01069 throw new Exception("\
01070 The options 1 to 5 don't work anymore. The reason is that now the code is\n\
01071 meshless and always assumes an orthonormal grid. This decision was made for\n\
01072 efficiency to save memory. So in this new version, all modules that used the\n\
01073 class Deal_II::Triangultion were deleted\n");
01074 }
01075 throw new Exception("Invalid option for DYNAMIC_MODULE\n");
01076 return NULL;
01077 }
01078 }
01079
01080 return m_pDynamicModule;
01081 }
01082
01083 void MainApp::printBeginSection(std::string str,std::ostream &out)
01084 {
01085 unsigned nSlashes = (70 - str.length())/2;
01086 out << "\n\n";
01087 for (unsigned i=0;i<nSlashes; i++)
01088 out << "=";
01089 printf(" %s ",str.c_str());
01090 for (unsigned i=0;i<nSlashes; i++)
01091 out << "=";
01092 out << "\n\n";
01093 }
01094
01095
01096
01097 void MainApp::printTriangulation()
01098 {
01099 std::string section="Mesh";
01100
01101 print(section,"Mesh: %g x %g x %g Elements\n",configFile.getDouble("NUM_ELEMS_X"),configFile.getDouble("NUM_ELEMS_Y"),configFile.getDouble("NUM_ELEMS_Z"));
01102 print(section,"Total: %d cells, %d faces\n",getMesh().numCells(),getMesh().numFaces());
01103 print(section,"Dimensions: %g x %g x %g\n",configFile.getDouble("DIM_X"),configFile.getDouble("DIM_Y"),configFile.getDouble("DIM_Z"));
01104 }
01105
01106
01107 void MainApp::textcolor(int attr, int fg, int bg)
01108 {
01109 char command[13];
01110
01111
01112 sprintf(command, "%c[%d;%d;%dm", 0x1B, attr, fg + 30, bg + 40);
01113 print("%s", command);
01114 }
01115
01116
01117
01124 void MainApp::getPBoundaryCondition(Function3D * &fDP, Function3D * &fFlux)
01125 {
01126 if (m_fDP)
01127 {
01128 fDP = m_fDP;
01129 fFlux = m_fFluxBC;
01130 return;
01131 }
01132
01133 std::string section="PBC";
01134
01135 print(section,"Boundary Condition for Pressure:\n");
01136 switch (configFile.getUnsigned("PRESSURE_BOUNDARY_CONDITIONS"))
01137 {
01138 case 1:
01139 {
01140 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01141 configFile.getDouble("DIM_Y"),
01142
01143
01144
01145
01146
01147
01148 configFile.getDouble("DIM_Z"),
01149 0,0,0,0,0,0);
01150
01151 fFlux =new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01152 configFile.getDouble("DIM_Y"),
01153 configFile.getDouble("DIM_Z"),
01154 INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY);
01155 print(section,"\tP(x) = 0 at boundary\n");
01156 break;
01157 }
01158 case 2:
01159 {
01160
01161 double dP = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01162 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01163 configFile.getDouble("DIM_Y"),
01164 configFile.getDouble("DIM_Z"),
01165 INFINITY,dP,INFINITY,
01166 INFINITY,INFINITY,INFINITY);
01167
01168 double dV = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01169 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01170 configFile.getDouble("DIM_Y"),
01171 configFile.getDouble("DIM_Z"),
01172 dV,INFINITY,0,0,0,0);
01173 print(section,"\tvel(x) = <%g,0> at left boundary\n\tvel(x) = 0 at top, bottom, front, back, boundaries\n\tp = %g at right boundary\n",dV,dP);
01174 break;
01175 }
01176 case 3:
01177 {
01178 double dPL = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01179 double dPR = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01180
01181 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01182 configFile.getDouble("DIM_Y"),
01183 configFile.getDouble("DIM_Z"),
01184 dPL,dPR,INFINITY,INFINITY,INFINITY,INFINITY);
01185 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01186 configFile.getDouble("DIM_Y"),
01187 configFile.getDouble("DIM_Z"),
01188 INFINITY,INFINITY,0.0,0.0,0.0,0.0);
01189 print(section,"\tP(x) = %g at left boundary,\n\tP(x) = %g at right boundary\n\tAll other boundaries are imperveous\n",dPL,dPR);
01190 break;
01191 }
01192 case 4:
01193 {
01194 Point3D DX = getMesh().getDX();
01195 double refP = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01196 fDP = new FCubeBoundaryCondition(DX(0),DX(1),DX(2),
01197 refP,INFINITY,INFINITY,
01198 INFINITY,INFINITY,INFINITY);
01199 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01200 configFile.getDouble("DIM_Y"),
01201 configFile.getDouble("DIM_Z"),
01202 0.0,0.0,0.0,0.0,0.0,0.0);
01203 fFlux = new FDomainComplement(fFlux,fDP);
01204 print(section,"\tVelocity equal to zero at boundaries\n");
01205 print(section,"\tReference pressure is %g setted in the bottom front left corner\n",refP);
01206 break;
01207 }
01208 case 5:
01209 {
01210 print(section,"Unidimensional pression boundary condition in X\t");
01211 double lp = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01212 double rp = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01213 print(section,"\nvel(x) = 0 at the top, bottom, front and back sides\n\t\
01214 p=%g, p=%g in the left and righ boundary\n",lp,rp);
01215
01216 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01217 configFile.getDouble("DIM_Y"),
01218 configFile.getDouble("DIM_Z"),
01219 lp,rp,INFINITY,INFINITY,INFINITY,INFINITY);
01220 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01221 configFile.getDouble("DIM_Y"),
01222 configFile.getDouble("DIM_Z"),
01223 INFINITY,INFINITY,0.0,0.0,0.0,0.0);
01224 break;
01225 }
01226
01227 case 6:
01228 {
01229 print(section,"Unidimensional pression boundary condition in Y\t");
01230 double lp = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01231 double rp = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01232 print(section,"vel(x) = 0 at the left, right , top and bottom sides\n\t\
01233 p=%g, p=%g in the front and back boundary\n",lp,rp);
01234
01235
01236 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01237 configFile.getDouble("DIM_Y"),
01238 configFile.getDouble("DIM_Z"),
01239 INFINITY,INFINITY,lp,rp,INFINITY,INFINITY);
01240 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01241 configFile.getDouble("DIM_Y"),
01242 configFile.getDouble("DIM_Z"),
01243 0.0,0.0,INFINITY,INFINITY,0.0,0.0);
01244 break;
01245 }
01246
01247
01248 case 7:
01249 {
01250 print(section,"Unidimensional pression boundary condition in Z\t");
01251 double bp = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01252 double up = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01253 print(section,"\nvel(x) = 0 at the left, right, front and back sides\n\t\
01254 p=%g, p=%g in the bottom and up boundary\n",bp,up);
01255
01256 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01257 configFile.getDouble("DIM_Y"),
01258 configFile.getDouble("DIM_Z"),
01259 INFINITY,INFINITY,INFINITY,INFINITY,bp,up);
01260 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01261 configFile.getDouble("DIM_Y"),
01262 configFile.getDouble("DIM_Z"),
01263 0.0,0.0,0.0,0.0,INFINITY,INFINITY);
01264 break;
01265 }
01266 case 8:
01267 {
01268 print(section,"Fivespot problem\n");
01269 double vI = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01270 double pO = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01271 double radius = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_3");
01272 print(section,"\tNormal Velocity at Injection Well: Vel(x)=%g\n",vI);
01273 print(section,"\tPressure at Production Well: Vel(x)=%g\n",pO);
01274 print(section,"\tWell radius: %g\n",radius);
01275 double height = getMesh().getQ()[2]-getMesh().getP()[2] + getMesh().getGridTolerance();
01276 Point3D C0 = getMesh().getP();
01277 C0[2]-=getMesh().getGridTolerance();
01278
01279
01280 Function3D *fAux0 = new FCylinderRegion(C0,radius,height,vI);
01281 Function3D *fAux1 = new FCylinderRegion(C0,radius,height,vI);
01282 C0[0]=getMesh().getQ()[0];
01283 C0[1]=getMesh().getQ()[1];
01284 fDP = new FCylinderRegion(C0,radius,height,pO);
01285 Function3D *fAux2 = new FCylinderRegion(C0,radius,height,pO);
01286
01287 fFlux=new FDomainUnion(fAux0,new FDomainComplement(new FPlane(0,0,0,0),new FDomainUnion(fAux1,fAux2),true),true);
01288
01289 break;
01290 }
01291 case 9:
01292 {
01293
01294 double dP = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_2");
01295 fDP = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01296 configFile.getDouble("DIM_Y"),
01297 configFile.getDouble("DIM_Z"),
01298 INFINITY,INFINITY,INFINITY,
01299 INFINITY,dP,INFINITY);
01300
01301 double dV = configFile.getDouble("PRESSURE_BOUNDARY_CONDITIONS_ARG_1");
01302 fFlux = new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01303 configFile.getDouble("DIM_Y"),
01304 configFile.getDouble("DIM_Z"),
01305 0,0,0,0,INFINITY,dV);
01306 print(section,"\tvel(x) = <-%g,0> at Upper boundaryy\n\tvel(x) = 0 at top, bottom, front, back, boundaries\n\tp = %g at bottom boundary\n",dV,dP);
01307 break;
01308 }
01309
01310
01311
01312 default:
01313 {
01314 throw new Exception("Invalid option for PRESSURE_BOUNDARY_CONDITION");
01315 }
01316 }
01317 Function3D *fPWellsBC,*fFluxWellsBC;
01318 getWellsPressureBC(fPWellsBC,fFluxWellsBC);
01319
01320 Function3D *f = new FDomainUnion(fFlux,fFluxWellsBC,true);
01321 fFlux = f;
01322
01323 f = new FDomainUnion(fDP,fPWellsBC,true);
01324 fDP=f;
01325
01326 m_fDP=fDP;
01327 m_fFluxBC=fFlux;
01328 }
01329
01330
01338 void MainApp::getUBoundaryCondition(Function3D * &fDU, Function3D * &fTensionBC)
01339 {
01340 std::string section="UBC";
01341
01342 print(section,"Boundary conditions for the displacement:\n");
01343 switch (configFile.getUnsigned("DISPLACEMENT_BOUNDARY_CONDITIONS"))
01344 {
01345 case 1:
01346 {
01347 fDU = new FCompoundVector( new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01348 configFile.getDouble("DIM_Y"),
01349 configFile.getDouble("DIM_Z"),
01350 0,0,INFINITY,INFINITY,INFINITY,INFINITY),
01351 new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01352 configFile.getDouble("DIM_Y"),
01353 configFile.getDouble("DIM_Z"),
01354 INFINITY,INFINITY,0,0,INFINITY,INFINITY),
01355 new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01356 configFile.getDouble("DIM_Y"),
01357 configFile.getDouble("DIM_Z"),
01358 INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,0),true);
01359
01360 double Weight = configFile.getDouble("DISPLACEMENT_BOUNDARY_CONDITIONS_ARG_1");
01361 fTensionBC = new FCompoundVector( new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01362 configFile.getDouble("DIM_Y"),
01363 configFile.getDouble("DIM_Z"),
01364 0,0,0,0,0,0),
01365 new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01366 configFile.getDouble("DIM_Y"),
01367 configFile.getDouble("DIM_Z"),
01368 0,0,0,0,0,0),
01369 new FCubeBoundaryCondition(configFile.getDouble("DIM_X"),
01370 configFile.getDouble("DIM_Y"),
01371 configFile.getDouble("DIM_Z"),
01372 0,0,0,0,0,-Weight),true);
01373 print(section,"\tThe domain is a vertical piston with a weight of %g at the top border\n",Weight);
01374 break;
01375 }
01376 default:
01377 {
01378 throw new Exception("Invalid option for PRESSURE_BOUNDARY_CONDITION");
01379 }
01380 }
01381 }
01382
01383
01387 Point3D MainApp::getSpatialDiscretizationSteps()
01388 {
01389 Point3D p;
01390 p[0] = configFile.getDouble("DIM_X")/configFile.getDouble("NUM_ELEMS_X");
01391 p[1] = configFile.getDouble("DIM_Y")/configFile.getDouble("NUM_ELEMS_Y");
01392 p[2] = configFile.getDouble("DIM_Z")/configFile.getDouble("NUM_ELEMS_Z");
01393
01394 return p;
01395 }
01396
01400 void MainApp::printMeshInfo()
01401 {
01402 if (configFile.getInt("DEBUG_PRINT_TRIANGULATION"))
01403 getMesh().print();
01404
01405 }
01406
01407
01408
01409
01410 Function3D* MainApp::getPorousMediaProperty(std::string key)
01411 {
01412 std::string section=key;
01413
01414 char str[1000];
01415 sprintf(str,"%s_ARG_1",key.c_str());
01416 std::string key_arg1=str;
01417
01418 sprintf(str,"%s_MEDIA",key.c_str());
01419 std::string key_media=str;
01420
01421 sprintf(str,"%s_FILE",key.c_str());
01422 std::string key_file=str;
01423
01424 sprintf(str,"%s_CV",key.c_str());
01425 std::string key_cv=str;
01426
01427 OrthoMesh &mesh = getMesh();
01428 print(section,"%s:\n",key.c_str());
01429 switch(configFile.getInt(key))
01430 {
01431 case 1:
01432 print(section,"\t%s = %g\n",key.c_str(),configFile.getDouble(key_arg1));
01433 return new FPlane(0.0,0.0,0.0,configFile.getDouble(key_arg1));
01434 break;
01435 case 2:
01436 {
01437 Point3D p0(0,0,0);
01438 Point3D p1(configFile.getDouble("DIM_X"),
01439 configFile.getDouble("DIM_Y"),
01440 configFile.getDouble("DIM_Z"));
01441 FChessBoard3D *f = new FChessBoard3D(p0,p1,mesh.getP(),mesh.getQ(),configFile.getString(key_file));
01442
01443 print(section,"\n\tChessBoard function\n\tfile %s\n\tDimensions: %d %d %d\n",
01444 configFile.getString(key_file).c_str(),
01445 f->getDimX(),f->getDimY(),f->getDimZ());
01446 return f;
01447 break;
01448 }
01449 case 3:
01450 {
01451 Point3D p0(0,0,0);
01452 Point3D p1(configFile.getDouble("DIM_X"),
01453 configFile.getDouble("DIM_Y"),
01454 configFile.getDouble("DIM_Z"));
01455 FChessBoard3D *f = new FChessBoard3D(p0,p1,mesh.getP(),mesh.getQ(),configFile.getDouble(key_media),
01456 configFile.getDouble(key_cv),
01457 configFile.getString(key_file));
01458
01459 print(section,"\n\tChessBoard with gaussian distribution function\n\tfile %s\n\tDimensions: %d %d %d\n\tMedia: %g\n\tCV=%g",
01460 configFile.getString(key_file).c_str(), f->getDimX(),f->getDimY(),f->getDimZ(),
01461 configFile.getDouble(key_media),configFile.getDouble(key_cv));
01462
01463 return f;
01464 break;
01465 }
01466
01467 default:
01468 {
01469 throw new Exception("Invalid option for key %s",key.c_str());
01470 return NULL;
01471 }
01472 }
01473 return NULL;
01474 }
01475
01476
01477
01478
01479
01480
01481
01482
01483
01489 FaceFluxFunction* MainApp::getFluxFunction()
01490 {
01491 std::string section="Flux";
01492 if (m_flux)
01493 {
01494 return m_flux;
01495 }
01496
01497 print(section,"Flux Function: ");
01498 switch(configFile.getInt("FLUX_FUNCTION"))
01499 {
01500 case 1:
01501 {
01502
01503 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
01504 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
01505
01506 print(section,"Flux function for biphasic flow without gravity term: F(u)=mobW Vdt\n");
01507
01508 m_flux= new FaceFluxWithoutGravity(getMesh(),dynamic_cast<Function1D&>(*fFracF),dynamic_cast<Function1D&>(*DfFracF));
01509 break;
01510 }
01511 case 2:
01512 {
01513 double v1 = configFile.getDouble("S1_VISCOSITY");
01514 double p1 = configFile.getDouble("S1_DENSITY");
01515 double v2 = configFile.getDouble("S2_VISCOSITY");
01516 double p2 = configFile.getDouble("S2_DENSITY");
01517 double g = configFile.getDouble("GRAVITY");
01518 print(section,"Bucley-Leverett Mobility with gravity term\n");
01519 print(section,"S1 Viscosity: %g\nS1 Density: %g\n",v1,p1);
01520 print(section,"S2 Viscosity: %g\nS2 Density: %g\n",v2,p2);
01521 print(section,"Gravity: %g\n",g);
01522
01523 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
01524 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
01525
01526
01527
01528 if (m_fK == NULL)
01529
01530 m_flux =new FaceFluxWithGravity(getMesh(),
01531 *dynamic_cast<Function1D*>(fFracF),
01532 *dynamic_cast<Function1D*>(DfFracF),
01533 *dynamic_cast<Function1D*>(fMobGrav),
01534 *dynamic_cast<Function1D*>(DfMobGrav),
01535 getPermeabilityAtCells(),-(p1-p2)*g);
01536 break;
01537 }
01538 case 3:
01539 print(section,"Burguers flux function\n");
01540 m_flux =new FaceFluxWithoutGravity(getMesh(),
01541 *(new FQuadratic(0.5,0.0,0.0)),
01542 *(new FLinear(1.0,0.0)));
01543 break;
01544 case 4:
01545 {
01546 print(section,"Independent system of equations\n\tBurguers for the first equation, Linear for the second equation:\n\t\
01547 F(S1,S2) = < S1*S1/2 * Vdt, S2* Vdt>\n\tWhere Vdt is the velocity of the dynamic module\n");
01548 m_flux= new FaceFluxWithIndependentEquations(getMesh(),
01549 new FQuadratic(0.5,0,0),
01550 new FLinear(1,0),
01551 new FLinear(1,0),
01552 new FLinear(0,1),true);
01553
01554 break;
01555 }
01556 case 5:
01557 {
01558
01559 double v1 = configFile.getDouble("S1_VISCOSITY");
01560 double v2 = configFile.getDouble("S2_VISCOSITY");
01561 double p1 = configFile.getDouble("S1_DENSITY");
01562 double p2 = configFile.getDouble("S2_DENSITY");
01563 double sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
01564 double sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
01565
01566 double g = configFile.getDouble("GRAVITY");
01567 print(section,"Flux designed for CO2 injection with dissolution using Bucley-Leverett Mobility and gravity term\n\
01568 F1(Sc,Cd) = fMob(Sc)*Vdt + K (1-C)(pw-pc) fMobGrav(Sc) g Grad(Z)\n\
01569 F2(Sc,Cd) = C*fMob(Sc)*Vdt - K C*(1-C)(pw-pc) fMobGrav(Sc) g Grad(Z)\n\
01570 \tWhere C = Cd/(1-Sc)\n\
01571 \tK = Permeability\n\
01572 \tpw = Water density\n\
01573 \tpc = CO2 density\n\
01574 \tCd = Is the fraction volume of the porous media occupied by the dissolved CO2\n\
01575 \tSc = Is the saturation of supercritic CO2\n\
01576 \tSc = Is the saturation of supercritic CO2\n\
01577 \tfMob(Sc) = Bucley Leverett relative mobility\n\
01578 \tfMobGrav(Sc) = The mobility for the gravity term\n\
01579 \tIn the program S1 represents supercritic CO2 while \nS2 is the dissolved volume of CO2 per porous media\n\
01580 CO2 Viscosity: %g\nCO2 Density (pc): %g\nCO2 Residual Saturation: %g\n\
01581 Water Viscosity: %g\nWater Density(pw): %g\nWater Residual Saturation: %g\n\
01582 Gravity: %g\n",v1,p1,sr1,v2,p2,sr2,g);
01583 m_flux= new FluxForCO2Inj(getMesh(),getPermeabilityAtCells(),v1,v2,sr1,sr2,p1,p2,g);
01584 break;
01585 }
01586 case 6:
01587 {
01588 double v1 = configFile.getDouble("S1_VISCOSITY");
01589 double v2 = configFile.getDouble("S2_VISCOSITY");
01590 print(section,"\
01591 Flux designed for Simple Black Oil Model of Oil and Gas where\n\
01592 Gas can dissolve into oil and oil does not evaporate.\n\
01593 Please see the article\n\
01594 \"John B Bell, Gregory R. Shubin, John Trangestein,\n\
01595 METHOD FOR REDUCING NUMERICAL DISPERSION IN TWO\n\
01596 PHASE BLACK OIL RESERVOIR SIMULATION, Journal of\n\
01597 Computational Physics 65, 71-106 1986\"\n\
01598 \n\
01599 The flux is a function of the total moles of each component\n \
01600 (i.e OIL mo, GAS mg using Bucley-Leverett Mobility.\n\
01601 The flux is given by the expression:\n\
01602 Fo(no,ng) = [m_o^o/volOil * fMobO(So) ]*Vdt\n\
01603 Fg(no,ng) = [m_g^o/volOil * fMobO(So) + m_g^g/VolGas * (1-fMobO(So))]*Vdt\n\
01604 Where\n\
01605 fMob = Bucley Leverret Mobility with\n\
01606 GAS Viscosity = Defined By Flash.\n\
01607 OIL Viscosity = Defined By Flash\n\
01608 So = Saturation of oil phase.\n\
01609 no,ng = Mass per pore volume of both components.\n\
01610 ",v1,v2);
01611 m_flux= new FaceFluxSimpleBlackOil(getMesh(),*getFlashCompositional());
01612 break;
01613 }
01614 case 7:
01615 {
01616 double v1 = configFile.getDouble("S1_VISCOSITY");
01617 double v2 = configFile.getDouble("S2_VISCOSITY");
01618 double sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
01619 double sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
01620 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
01621 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
01622
01623
01624 double grav = configFile.getDouble("GRAVITY");
01625 print(section,"\
01626 Flux designed for CO2 injection in Brine for compositional model.\n\
01627 The flux is a function of the total moles of each component\n\
01628 (i.e Water,CO2, SALT) using Bucley-Leverett Mobility.\n\
01629 The flux is given by the expression:\n\
01630 F1(mw,mc,ms) = [mw_aq/volAqueous * VelAQ + mw_gas/VolGas * VelG\n\
01631 F2(mw,mc,ms) = [mc_aq/volAqueous * VelAQ + mc_gas/VolGas * VelG\n\
01632 F3(mw,mc,ms) = [ms_aq/volAqueous * VelAQ\n\
01633 Where\n\
01634 VelAQ = mobW(Sw)*Vdt + K*mobW*(1-mobW)*mobT(densityG - densityAQ)*grav \n\
01635 VelO = Vdt-VelAQ\n\
01636 K = Permeability Feld\n\
01637 fMob = Bucley Leverret Mobility with\n\
01638 Water Viscosity = %g.\n\
01639 CO2 Viscosity = %g.\n\
01640 Water Residual Saturation = %g\n\
01641 CO2 Residual Saturation = %g\n\
01642 Gravity = %g.\n\
01643 Sw = Saturation of water phase.\n\
01644 mw,mc,mw = Moles per pore volume of each component.\n"
01645 ,v1,v2,sr1,sr2,grav);
01646 m_flux= new FaceFluxCO2Brine(getMesh(),*getFlashCompositional(),getPermeabilityAtCells(),*fFracF,*DfFracF,*fMobGrav,*DfMobGrav,grav);
01647 break;
01648 }
01649 case 8:
01650 {
01651 print(section,"Flux Mass Simple Black Oil Model\n");
01652 m_flux= new FaceFluxSimpleBlackOilMass(getMesh(),*getFlashCompositional());
01653 break;
01654 }
01655
01656 default:
01657 {
01658 throw new Exception("Invalid Option FLUX_FUNCTION=%d\n",
01659 configFile.getInt("FLUX_FUNCTION"));
01660 return NULL;
01661 }
01662 }
01663 return m_flux;
01664 }
01665
01666 VecWellInfo MainApp::getWells()
01667 {
01668 static bool firstTime=true;
01669
01670 if (!firstTime)
01671 return m_wells;
01672 firstTime=false;
01673 VecWellInfo &wells=m_wells;
01674 unsigned i=1;
01675 char keyName[500];
01676 std::string section="Wells";
01677 while (1)
01678 {
01679 sprintf(keyName,"WELL_%d_PX",i);
01680 if (configFile.isDefined(keyName))
01681 {
01682
01683 Point3D P;
01684 sprintf(keyName,"WELL_%d_PX",i);
01685 P(0)=configFile.getDouble(keyName);
01686 sprintf(keyName,"WELL_%d_PY",i);
01687 P(1)=configFile.getDouble(keyName);
01688 sprintf(keyName,"WELL_%d_PZ",i);
01689 P(2)=configFile.getDouble(keyName);
01690
01691 Point3D Q;
01692 sprintf(keyName,"WELL_%d_QX",i);
01693 Q(0)=configFile.getDouble(keyName);
01694 sprintf(keyName,"WELL_%d_QY",i);
01695 Q(1)=configFile.getDouble(keyName);
01696 sprintf(keyName,"WELL_%d_QZ",i);
01697 Q(2)=configFile.getDouble(keyName);
01698
01699 sprintf(keyName,"WELL_%d_INJECTION_RATE",i);
01700 if (configFile.isDefined(keyName))
01701 {
01702 double value = configFile.getDouble(keyName);
01703 wells.push_back(WellInfo(P,Q,value,WellInfo::FLUX));
01704 }
01705 else
01706 {
01707 sprintf(keyName,"WELL_%d_PRESSURE",i);
01708 if (configFile.isDefined(keyName))
01709 {
01710 double value = configFile.getDouble(keyName);
01711 wells.push_back(WellInfo(P,Q,value,WellInfo::PRESSURE));
01712 }
01713 else
01714 {
01715 sprintf(keyName,"WELL_%d_SOURCE_INJECTION",i);
01716 if (configFile.isDefined(keyName))
01717 {
01718 double value = configFile.getDouble(keyName);
01719 wells.push_back(WellInfo(P,Q,value,WellInfo::SOURCE_INJ_RATE));
01720 }
01721 else
01722 {
01723 sprintf(keyName,"WELL_%d_SOURCE_FIXED_PRESSURE",i);
01724 if (configFile.isDefined(keyName))
01725 {
01726 double value = configFile.getDouble(keyName);
01727 wells.push_back(WellInfo(P,Q,value,WellInfo::SOURCE_FIXED_PRESSURE));
01728 }
01729 else
01730 throw new Exception("Expect one of the keys WELL_%d_SOURCE_INJECTION WELL_%d_SOURCE_FIXED_PRESSURE WELL_%d_PRESSURE WELL_%d_INJECTION_RATE",i,i,i,i);
01731 }
01732 unsigned j=1;
01733 sprintf(keyName,"WELL_%d_S%d_BC",i,j++);
01734 std::vector<double> v1;
01735 while(configFile.isDefined(keyName))
01736 {
01737 v1.push_back(configFile.getDouble(keyName));
01738 sprintf(keyName,"WELL_%d_S%d_BC",i,j++);
01739 }
01740 if (v1.size() < getFluxFunction()->numComponents())
01741 throw new Exception("Error, Number of transport boundary conditions for the well %d is %d but transport has %d variables",i,v1.size(),getFluxFunction()->numComponents());
01742 wells.back().setTransportBC(v1);
01743 }
01744 }
01745 }
01746 else
01747 break;
01748 i++;
01749 }
01750
01751 Point3D DX = getDX();
01752 double tol = NumericMethods::min(DX[0],DX[1],DX[2])/10.0;
01753
01754 if (!wells.size())
01755 print(section, "No Wells were specified\n");
01756 for (unsigned i=0;i<wells.size();i++)
01757 {
01758 wells[i].adjustBoundaryWithGrid(DX);
01759
01760
01761
01762 OrthoMesh &mesh = getMesh();
01763 if (NumericMethods::a_equal(wells[i].P[0],mesh.getP()[0],tol) ||
01764 NumericMethods::a_equal(wells[i].P[1],mesh.getP()[1],tol) ||
01765 NumericMethods::a_equal(wells[i].P[2],mesh.getP()[2],tol) ||
01766 NumericMethods::a_equal(wells[i].Q[0],mesh.getQ()[0],tol) ||
01767 NumericMethods::a_equal(wells[i].Q[1],mesh.getQ()[1],tol) ||
01768 NumericMethods::a_equal(wells[i].Q[2],mesh.getQ()[2],tol))
01769 throw new Exception("Well %d intersects the boundary of the domain or subdomain");
01770
01771 if (wells[i].getBCType() == WellInfo::FLUX)
01772 {
01773 print(section,"Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tInjection Rate: %g\n",i+1,
01774 wells[i].P[0],wells[i].P[1],wells[i].P[2],
01775 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
01776 wells[i].volume(),wells[i].getInjRate());
01777 }
01778 else if (wells[i].getBCType() == WellInfo::PRESSURE)
01779 {
01780 print(section,"Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tPressure: %g\n\t",i+1,
01781 wells[i].P[0],wells[i].P[1],wells[i].P[2],
01782 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
01783 wells[i].volume(),wells[i].getPressureBC());
01784 }
01785 else if (wells[i].getBCType() == WellInfo::SOURCE_INJ_RATE )
01786 {
01787 print(section,"Source Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tTotal Injection: %g\n\t",i+1,
01788 wells[i].P[0],wells[i].P[1],wells[i].P[2],
01789 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
01790 wells[i].volume(),wells[i].getSourceInjRate());
01791 const VecDouble &values=wells[i].getTransportBC();
01792 for (unsigned j=0;j<values.size();j++)
01793 {
01794 print(section,"S%d = %g\n\t",j,values(j));
01795 }
01796 }
01797 else if (wells[i].getBCType() == WellInfo::SOURCE_FIXED_PRESSURE )
01798 {
01799 print(section,"Source Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tFixed Pressure: %g\n\t",i+1,
01800 wells[i].P[0],wells[i].P[1],wells[i].P[2],
01801 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
01802 wells[i].volume(),wells[i].getSourceFixedPressure());
01803 const VecDouble &values=wells[i].getTransportBC();
01804 for (unsigned j=0;j<values.size();j++)
01805 {
01806 print(section,"S%d = %g\n\t",j,values(j));
01807 }
01808 }
01809
01810
01811 }
01812 return wells;
01813 }
01814
01815
01816 VecWellInfo MainApp::getSourceWells()
01817 {
01818 VecWellInfo wells=getWells();
01819 VecWellInfo sourceWells;
01820 for (unsigned i=0;i<wells.size();i++)
01821 {
01822 if (wells[i].getBCType() == WellInfo::SOURCE_INJ_RATE)
01823 {
01824 sourceWells.push_back(wells[i]);
01825 }
01826 }
01827 return sourceWells;
01828
01829 }
01830
01831
01832
01833 Point3D MainApp::getDomainDimensions(){
01834 return Point3D(configFile.getDouble("DIM_X"),
01835 configFile.getDouble("DIM_Y"),
01836 configFile.getDouble("DIM_Z"));
01837 }
01838
01839
01840
01841
01842 VecDouble& MainApp::getPorosityAtCells()
01843 {
01844 if (m_cPor.size() == 0)
01845 {
01846 m_cPor.reinit(getMesh().numCells());
01847 m_fPor = getPorousMediaProperty("POROSITY");
01848 getMesh().setCentralValuesFromFunction(*m_fPor,m_cPor);
01849 }
01850 return m_cPor;
01851 }
01852
01853 VecDouble& MainApp::mb_getPorosityAtCells()
01854 {
01855 if (m_mbcPor.size() == 0)
01856 {
01857 m_mbcPor.reinit(getMesh().numCells());
01858 m_mbfPor = getPorousMediaProperty("MB_POROSITY");
01859 getMesh().setCentralValuesFromFunction(*m_mbfPor,m_mbcPor);
01860 }
01861 return m_mbcPor;
01862 }
01863
01864
01865
01866 VecDouble& MainApp::getPermeabilityAtCells()
01867 {
01868 if (m_cK.size() == 0)
01869 {
01870 m_cK.reinit(getMesh().numCells());
01871 Function3D *f = getPorousMediaProperty("PERMEABILITY");
01872 getMesh().setCentralValuesFromFunction(*f,m_cK);
01873 delete f;
01874
01875
01876
01877 VecDouble vV(getMesh().numVertices());
01878 getMesh().projectCentralValuesAtVertices(m_cK,vV);
01879 HDF5OrthoWriter::getHDF5OrthoWriter().writeScalarField(vV,"K");
01880 }
01881 return m_cK;
01882
01883 }
01884
01885 VecDouble& MainApp::mb_getPermeabilityAtCells()
01886 {
01887 if (m_mbcK.size() == 0)
01888 {
01889 m_mbcK.reinit(getMesh().numCells());
01890 Function3D *fmb = getPorousMediaProperty("MB_PERMEABILITY");
01891 getMesh().setCentralValuesFromFunction(*fmb,m_mbcK);
01892 delete fmb;
01893
01894 }
01895 return m_mbcK;
01896
01897 }
01898
01899
01900
01901
01902
01903 ConstVelocityModule* MainApp::getConstVelocityDynamicModule()
01904 {
01905 std::string section="ConstVelocity";
01906 Function3D *fP = getPressureIC();
01907 switch (configFile.getInt("CONST_VELOCITY_FUNCTION"))
01908 {
01909 case 1:
01910 print(section,"Vel(x) = <1,0,0>\n");
01911 this->m_pVelFunction = new ConstVectorFunction(1,0,0);
01912 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01913 case 2:
01914 print(section,"Vel(x) = <-1,0,0>\n");
01915 this->m_pVelFunction = new ConstVectorFunction(-1,0,0);
01916 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01917 case 3:
01918 print(section,"Vel(x) = <0,1,0>\n");
01919 this->m_pVelFunction = new ConstVectorFunction(0,1,0);
01920 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01921 case 4:
01922 print(section,"Vel(x) = <0,-1,0>\n");
01923 this->m_pVelFunction = new ConstVectorFunction(0,-1,0);
01924 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01925 case 5:
01926 print(section,"Vel(x) = <0,0,1>\n");
01927 this->m_pVelFunction = new ConstVectorFunction(0,0,1);
01928 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01929 case 6:
01930 print(section,"Vel(x) = <0,0,-1>\n");
01931 this->m_pVelFunction = new ConstVectorFunction(0,0,-1);
01932 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01933 case 7:
01934 print(section,"Vel(x) = <1,0,0> p/ y<YDIM/2\nVel(x) = <0.2,0,0> y > YDIM/2\n");
01935 this->m_pVelFunction = new TwoTracksVectorFunction(configFile.getDouble("DIM_Y")/2.0,1,0.2,X,Y);
01936 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01937 case 8:
01938 print(section,"Vel(x) = <0,0,0>\n");
01939 this->m_pVelFunction = new ConstVectorFunction(0,0,0);
01940 return new ConstVelocityModule(getMesh(),*m_pVelFunction,*fP);
01941 default:
01942 throw new Exception("Invalid option for CONST_VELOCITY_FUNCTION");
01943 return NULL;
01944 }
01945 }
01946
01947
01948
01949
01950
01951
01952
01953
01954 unsigned MainApp::hasSystemOfTransportEquations()
01955 {
01956 enablePrint(false);
01957 FaceFluxFunction *pFlux = getFluxFunction();
01958 enablePrint(true);
01959
01960 if (pFlux == NULL)
01961 return true;
01962 else
01963 {
01964 delete pFlux;
01965 return false;
01966 }
01967
01968 }
01969
01970
01971
01972
01973
01974
01975
01976
01977
01981 Point3D MainApp::getP0()
01982 {
01983 return Point3D(0.0,0.0,0.0);
01984 }
01985
01989 Point3D MainApp::getP1()
01990 {
01991 return Point3D(configFile.getDouble("DIM_X"),
01992 configFile.getDouble("DIM_Y"),
01993 configFile.getDouble("DIM_Z"));
01994 }
01995
01996
02000 Point3D MainApp::getDX()
02001 {
02002 Point3D p;
02003 p[0] = configFile.getDouble("DIM_X")/configFile.getDouble("NUM_ELEMS_X");
02004 p[1] = configFile.getDouble("DIM_Y")/configFile.getDouble("NUM_ELEMS_Y");
02005 p[2] = configFile.getDouble("DIM_Z")/configFile.getDouble("NUM_ELEMS_Z");
02006
02007 return p;
02008 }
02009
02010
02011 void MainApp::getWellsPressureBC(Function3D *&fP,Function3D *&fFlux)
02012 {
02013 VecWellInfo wells = getWells();
02014 VecWellInfo fluxWells;
02015 VecWellInfo pWells;
02016
02017 std::vector<double> fluxValues;
02018 std::vector<double> pValues;
02019 Point3D DX = getDX();
02020 DX/=4.0;
02021
02022 std::string section="Wells";
02023
02024 for (unsigned i=0;i<wells.size();i++)
02025 {
02026 if (wells[i].getBCType() == WellInfo::FLUX)
02027 {
02028 print(section,"Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tInjection Rate: %g\n\t",i,
02029 wells[i].P[0],wells[i].P[1],wells[i].P[2],
02030 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
02031 wells[i].volume(),wells[i].getInjRate());
02032 wells[i].P-=DX;
02033 wells[i].Q+=DX;
02034 fluxWells.push_back(wells[i]);
02035 fluxValues.push_back(wells[i].getInjRate()/wells[i].area());
02036 print(section,"Velocity at faces: %g\n",wells[i].getInjRate()/wells[i].area());
02037
02038 }
02039 else
02040 {
02041 print(section,"Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\tPressure: %g\n\t",i,
02042 wells[i].P[0],wells[i].P[1],wells[i].P[2],
02043 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
02044 wells[i].volume(),wells[i].getPressureBC());
02045 wells[i].P-=DX;
02046 wells[i].Q+=DX;
02047 pValues.push_back(wells[i].getPressureBC());
02048 pWells.push_back(wells[i]);
02049
02050 }
02051 }
02052 VecDouble vFluxValues;
02053 VecDouble vPValues;
02054 NumericMethods::STLToVecDouble(vFluxValues,fluxValues);
02055 NumericMethods::STLToVecDouble(vPValues,pValues);
02056 fP= new FWellCondition(pWells,vPValues);
02057 fFlux=new FWellCondition(fluxWells,vFluxValues);
02058 }
02059
02060 Function3D& MainApp::getWellsSourceTerm()
02061 {
02062 if (m_fSource)
02063 return *m_fSource;
02064
02065 VecWellInfo wells = getWells();
02066 VecWellInfo sourceWells;
02067 std::vector<double> values;
02068 Point3D DX = getDX();
02069 DX/=4.0;
02070
02071 std::string section="SourceWells";
02072
02073 for (unsigned i=0;i<wells.size();i++)
02074 {
02075 if (wells[i].getBCType() == WellInfo::SOURCE_INJ_RATE)
02076 {
02077 print(section,"Source Well %d => <%g,%g,%g> - <%g,%g,%g>\n\tvolume = %g\n\t\n\tarea: %g\n\tTotal Flux Rate: %g\n\t",i,
02078 wells[i].P[0],wells[i].P[1],wells[i].P[2],
02079 wells[i].Q[0],wells[i].Q[1],wells[i].Q[2],
02080 wells[i].volume(),wells[i].area(),wells[i].getSourceInjRate());
02081 wells[i].P-=DX;
02082 wells[i].Q+=DX;
02083 sourceWells.push_back(wells[i]);
02084 }
02085 }
02086 m_fSource=new FWellsSource(sourceWells);
02087 return *m_fSource;
02088 }
02089
02090
02091 Function3D* MainApp::getWellTransportBC(std::string varName)
02092 {
02093
02094 std::string section="WellTransportBC";
02095
02096
02097 VecWellInfo wells = getWells();
02098 VecDouble values(wells.size());
02099 Point3D DX = getDX();
02100 DX/=4.0;
02101
02102
02103 for (unsigned i=0;i<wells.size();i++)
02104 {
02105 char key[100];
02106 sprintf(key,"WELL_%d_%s_BC",i+1,varName.c_str());
02107 double dd = configFile.getDouble(key);
02108 print(section,"Well %d => %s = %g\n",i+1,varName.c_str(),dd);
02109 values(i)=dd;
02110
02111 wells[i].P-=DX;
02112 wells[i].Q+=DX;
02113
02114
02115 }
02116 return new FWellCondition(wells,values);
02117
02118 }
02119
02120
02121 void MainApp::print(std::string section,const char *format,...)
02122 {
02123 char strOut[5000];
02124
02125 if (m_bPrint)
02126 {
02127 va_list vl;
02128 va_start(vl,format);
02129 vsprintf(strOut,format,vl);
02130 va_end(vl);
02131 log[section]+=strOut;
02132
02133
02134 }
02135
02136
02137 }
02138
02139
02140 std::vector<unsigned> MainApp::getNElements()
02141 {
02142 std::vector<unsigned> v(3);
02143 v[0]=(unsigned) configFile.getInt("NUM_ELEMS_X");
02144 v[1]=(unsigned) configFile.getInt("NUM_ELEMS_Y");
02145 v[2]=(unsigned) configFile.getInt("NUM_ELEMS_Z");
02146 return v;
02147 }
02148
02149
02150
02151 class CoarseOp : public HDF5Operator
02152 {
02153 private:
02154 HDF5OrthoReader &orthoReader;
02155 HDF5OrthoWriter &orthoWriter;
02156 unsigned nX,nY,nZ;
02157 public:
02158 CoarseOp(HDF5OrthoReader &orthoReader,HDF5OrthoWriter &orthoWriter,unsigned nX,unsigned nY,unsigned nZ)
02159 :orthoReader(orthoReader),orthoWriter(orthoWriter),nX(nX),nY(nY),nZ(nZ)
02160 {
02161
02162 }
02163
02164 virtual void processDataSet(hid_t dataset)
02165 {
02166 char dataSetName[1000];
02167
02168 H5Iget_name(dataset,dataSetName,1000);
02169
02170 ArrayString strLst=StringFunctions::tokenize(dataSetName,"/");
02171 if (strLst.size() == 0)
02172 return;
02173 if (strLst[0] == "Triangulations")
02174 {
02175 printf("Ignoring dataset %s\n",dataSetName);
02176 return;
02177 }
02178
02179
02180 printf("Coping DataSet %s\n",dataSetName);
02181 static VecDouble v;
02182
02183
02184 orthoReader.readScalars(v,dataSetName);
02185
02186
02187 std::string triaName = orthoReader.readAtt(dataSetName,"Triangulation");
02188 if (triaName.empty())
02189 throw new Exception("Error processing DataSet %s, Attribute \"Triangulation\" is missing\n",dataSetName);
02190 if (!orthoWriter.existTriaInformation(triaName))
02191 {
02192
02193 printf("===============Reading Mesh %s...",triaName.c_str());
02194 OrthoMesh &mesh = orthoReader.getMesh(triaName);
02195 printf("Ok\n");
02196 printf("===============Writing Mesh %s...",triaName.c_str());
02197 orthoWriter.writeCoarseMesh(triaName,mesh,nX,nY,nZ);
02198 printf("Ok\n");
02199
02200 }
02201
02202 printf("===============Writing DataSet %s...",dataSetName);
02203 orthoWriter.writeDataSet(v,dataSetName,triaName);
02204 printf("Ok\n");
02205
02206
02207 }
02208 virtual void processAtt(hid_t att)
02209 {
02210 char str[1000];
02211 char str2[1000];
02212
02213 H5Iget_name(att,str,1000);
02214 H5Aget_name(att,1000,str2);
02215
02216 ArrayString strLst=StringFunctions::tokenize(str,"/");
02217 if (strLst.size() == 0)
02218 return;
02219 if (strLst[0] == "Triangulations")
02220 {
02221 printf("Ignoring Attribute %s/%s\n",str,str2);
02222 return;
02223 }
02224
02225
02226
02227 string attValue=orthoReader.readAtt(str,str2);
02228 printf("Coping attribute %s/%s = %s\n",str,str2,attValue.c_str());
02229 orthoWriter.setAtt(str,str2,attValue);
02230 }
02231
02232 virtual void processGroup(hid_t grp)
02233 {
02234 char str[1000];
02235 H5Iget_name(grp,str,1000);
02236 ArrayString strLst=StringFunctions::tokenize(str,"/");
02237 if (strLst.size() == 0)
02238 return;
02239 if (strLst[0] != "Triangulations")
02240 {
02241 printf("Creating Group %s\n",str);
02242 orthoWriter.createGroup(str);
02243 }
02244 }
02245
02246
02247
02248 };
02249
02250
02251
02252 void MainApp::coarseMesh(unsigned nX,unsigned nY,unsigned nZ,string fileIn,string fileOut)
02253 {
02254 try{
02255
02256 printBeginSection("Data File Coarsing option");
02257 printf("Blocks in X: %d\n",nX);
02258 printf("Blocks in Y: %d\n",nY);
02259 printf("Blocks in Z: %d\n",nZ);
02260 printf("Input File: %s\n",fileIn.c_str());
02261 printf("Output File: %s\n",fileOut.c_str());
02262
02263
02264 HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
02265 hdf5.setOutputFile(fileOut);
02266 HDF5OrthoReader hdf5In;
02267 hdf5In.open(fileIn);
02268 CoarseOp op(hdf5In,HDF5OrthoWriter::getHDF5OrthoWriter(),nX,nY,nZ);
02269 hid_t root = H5Gopen1(hdf5In.getFile(),"/");
02270 hdf5In.iterateTree(root,op);
02271 H5Gclose(root);
02272 }
02273 catch(Exception *e)
02274 {
02275 printf("Error");
02276 printBeginSection("ERROR MESSAGE");
02277 printf("%s\n\nQuitting.....\n",e->getError().c_str());
02278 return;
02279 }
02280 }
02281
02282
02283
02284
02285
02286
02287 double MainApp::getYoung()
02288 {
02289 double result = configFile.getDouble("YOUNG_COEFFICIENT");
02290 print("YOUNG","YOUNG COEFFICIENT = %g\n",result);
02291 return result;
02292 }
02293
02294 double MainApp::getPoisson()
02295 {
02296 double result = configFile.getDouble("POISSON_COEFFICIENT");
02297 print("POISSON","POISSON COEFFICIENT = %g\n",result);
02298 return result;
02299 }
02300
02301
02302
02303 double MainApp::getDynamicTimeStep()
02304 {
02305 double nIt = configFile.getDouble("N_DYN_ITERATIONS");
02306 double time = configFile.getDouble("END_TIME");
02307 return time/nIt;
02308 }
02309
02310
02311 LinearSolver* MainApp::getLinearSolver()
02312 {
02313 static char str[500];
02314 if (!m_pSolver)
02315 {
02316 unsigned solverId = configFile.getInt("LINEAR_SOLVER");
02317 switch (solverId)
02318 {
02319 case 1:
02320 sprintf(str,"Solver: UMFPACK (Direct Solver for Sparse Matrix)\n");
02321 m_pSolver= new UMFPACKSolver();
02322 break;
02323 case 2:
02324 {
02325 double tol = configFile.getDouble("SOLVER_TOLERANCE");
02326 unsigned nIt = configFile.getInt("SOLVER_NUM_ITERATIONS");
02327 sprintf(str,"Solver: CG with Algebraic Multigrid Pre-conditioning (AGM)\nfor Sparse Symetric Positive Definite Matrices\n\tTolerance: %g\n\tMax Num Iterations: %d\n",tol,nIt);
02328 m_pSolver= new SolverGPU_AGM(tol,nIt);
02329 break;
02330 }
02331 case 3:
02332 {
02333 double tol = configFile.getDouble("SOLVER_TOLERANCE");
02334 unsigned nIt = configFile.getInt("SOLVER_NUM_ITERATIONS");
02335 sprintf(str,"Solver: CG with SSOR preconditioner for Symetric Definite Matrices\n\tTolerance: %g\n\tMax Num Iterations: %d\n",tol,nIt);
02336 unsigned debugLevel = configFile.getInt("DYNAMIC_DEBUG_LEVEL");
02337 m_pSolver= new CGSolver(nIt,tol,debugLevel);
02338 break;
02339 }
02340
02341 default:
02342 throw new Exception("Invalid Option (%d) for the Linear Solver\n",configFile.getInt("LINEAR_SOLVER"));
02343 return NULL;
02344 }
02345 }
02346 print("LINEAR_SOLVER",str);
02347 return m_pSolver;
02348 }
02349
02350
02360 void MainApp::adjustInitialConditionForCompressibleModel(DynamicBase &dyn,ConservativeMethodForSystem &trans, FlashCompositional &flash)
02361 {
02362 flash.setTransport(trans);
02363 flash.setDynamic(dyn);
02364 unsigned nCells=getMesh().numCells();
02365 unsigned nComps=flash.numComponents();
02366 VecDouble data(m_flux->numComponents());
02367 VecDouble phasesVol(flash.numPhases());
02368 FlashData flashData(flash.numPhases(),flash.numComponents(),NULL);
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380 flash.execute();
02381 for (unsigned cell=0;cell<nCells;cell++)
02382 {
02383 flash.flash(cell,flashData);
02384 flash.getPhasesVolume(dyn.getPressureAtCells()(cell),flashData,phasesVol);
02385 double factor=1.0/NumericMethods::vectorSum(phasesVol);
02386
02387
02388 for (unsigned i=0;i<nComps;i++)
02389 trans.getSolutionAtCells()(cell,i)*=factor;
02390 }
02391 flash.execute();
02392
02393 Function3D *fDP,*fFlux;
02394 getPBoundaryCondition(fDP,fFlux);
02395 trans.updateFixedPressureBC(fDP,flash);
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408 }
02409
02410
02411
02416 Function3D* MainApp::getTransportForSystemBC(unsigned nCmps)
02417 {
02418 std::vector<Function3D*> fV;
02419 for (unsigned i=0;i<nCmps;i++)
02420 {
02421 char str[10];
02422 sprintf(str,"S%d",i+1);
02423 fV.push_back(getTransportBC(str));
02424 }
02425 return new FCompoundVector(fV);
02426 }
02427
02432 Function3D* MainApp::getTransportForSystemIC(unsigned nCmps)
02433 {
02434 std::vector<Function3D*> fV;
02435 for (unsigned i=0;i<nCmps;i++)
02436 {
02437 char str[10];
02438 sprintf(str,"S%d",i+1);
02439 fV.push_back(getTransportIC(str));
02440 }
02441 return new FCompoundVector(fV);
02442 }
02443
02444
02445 FixedValueCondition& MainApp::getTransportFixedConditions()
02446 {
02447 static bool firstTime=true;
02448 if (firstTime){
02449
02450 std::string section="FixedSi";
02451 unsigned i=1;
02452 char key_name[200]="S1_FIXED_1_X";
02453 Point3D p;
02454 OrthoMesh &mesh = getMesh();
02455 if (configFile.isDefined(key_name))
02456 print(section,"Fixed Boundary Conditions:\n");
02457
02458 while (configFile.isDefined(key_name))
02459 {
02460
02461 sprintf(key_name,"S1_FIXED_%d_X",i);
02462 p[0]=configFile.getDouble(key_name);
02463
02464 sprintf(key_name,"S1_FIXED_%d_Y",i);
02465 p[1]=configFile.getDouble(key_name);
02466
02467 sprintf(key_name,"S1_FIXED_%d_Z",i);
02468 p[2]=configFile.getDouble(key_name);
02469
02470 sprintf(key_name,"S1_FIXED_%d_VALUE",i);
02471 double value = configFile.getDouble(key_name);
02472
02473 m_fixedC.addFixedCondition(mesh,p,value);
02474
02475 print(section,"\tS1(%g,%g,%g) = %g\n",p[0],p[1],p[2],value);
02476
02477 i++;
02478 sprintf(key_name,"S1_FIXED_%d_X",i);
02479 }
02480
02481
02482
02483 VecWellInfo wells = getSourceWells();
02484
02485
02486
02487
02488
02489 m_fixedC. addTransportFixedCondition(mesh,wells);
02490 firstTime=false;
02491 }
02492 return m_fixedC;
02493 }
02494
02495 FixedValueCondition& MainApp::getPressureFixedConditions()
02496 {
02497 static bool firstTime=true;
02498 if (firstTime){
02499 firstTime=false;
02500 VecWellInfo wells = getSourceWells();
02501 OrthoMesh::Cell_It cell = getMesh().begin_cell();
02502 OrthoMesh::Cell_It endc = getMesh().end_cell();
02503 for(;cell!=endc;cell++)
02504 {
02505
02506 Point3D p;
02507 cell->barycenter(p);
02508 for (unsigned i=0;i<wells.size();i++)
02509 {
02510
02511 if (wells[i].getBCType() == WellInfo::SOURCE_FIXED_PRESSURE && wells[i].isPointInWell(p))
02512 {
02513 m_fixedP.addFixedCondition(cell->index(),wells[i].getSourceFixedPressure());
02514 }
02515 }
02516 }
02517 }
02518
02519
02520 return m_fixedP;
02521 }
02522
02523 std::string MainApp::getHDF5FileName()
02524 {
02525 if (configFile.isDefined("HDF5_OUTPUT"))
02526 return configFile.getString("HDF5_OUTPUT");
02527 else
02528 {
02529 string file=configFile.getCurrentFile();
02530 size_t start=file.find_last_of('/');
02531 if (start== string::npos)
02532 start = 0;
02533 else
02534 start++;
02535
02536 size_t end=file.find_last_of('.');
02537 if (end == string::npos)
02538 end=file.size();
02539
02540
02541 string resp=file.substr(start,end-start);
02542 resp+=m_suffix+".hdf";
02543 return resp;
02544 }
02545 }
02546 using namespace StringFunctions;
02547 std::string MainApp::getMonitorWellsOutputFileName()
02548 {
02549 string file=configFile.getCurrentFile();
02550 std::string resp=get_prefix_file(get_file_name(configFile.getCurrentFile()));
02551 return "mw_" + resp + m_suffix+".dat";
02552 }
02553
02554
02555
02556
02557
02558
02559 void MainApp::abortEvent()
02560 {
02561
02562 if (m_pDynamicModule)
02563 m_pDynamicModule->printOutput();
02564 if (m_pTransportModule)
02565 m_pTransportModule->printOutput();
02566 if (m_pFlash)
02567 m_pFlash->printOutput();
02568 }
02569
02570
02571 BlockMatrixModule& MainApp::getBlockMatrix()
02572 {
02573
02574 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
02575 mb_getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
02576
02577 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
02578 mb_getPcFunction(&fPc,&dfPc,&fPcInv);
02579
02580 VecDouble pcIC(getMesh().numCells());
02581
02582
02583
02584 Function3D *f = getTransportIC("MB_S1");
02585 this->getMesh().setCentralValuesFromFunction(*f,pcIC);
02586 (dynamic_cast<Function1D&>(*fPc)).evaluate(pcIC);
02587
02588
02589
02590 return *(new BlockMatrixModule(this->getMB_Mesh() ,getMesh().numCells(), pcIC,mb_getPorosityAtCells(), mb_getPermeabilityAtCells(),
02591 dynamic_cast<Function1D&>(*dfPc),
02592 dynamic_cast<Function1D&>(*fPcInv),
02593 dynamic_cast<Function1D&>(*fMobT),
02594 dynamic_cast<Function1D&>(*fMob)));
02595
02596 }
02597
02598
02599 DiffusiveStep* MainApp::getDiffusiveStep()
02600 {
02601 std::string section="DiffusiveStep";
02602 if (m_pDiff)
02603 return m_pDiff;
02604
02605 switch(configFile.getInt("DIFFUSIVE_STEP",0))
02606 {
02607 case 0:
02608 {
02609
02610 if (configFile.getInt("DYNAMIC_MODULE")==7)
02611 {
02612 throw new Exception("Diffusive Step must be setted for dynamic module 7\n");
02613 }
02614 m_pDiff=NULL;
02615 break;
02616 }
02617 case 1:
02618 {
02619 print(section,"Diffusive Step\n");
02620 print(section,"\tMixed Hybrid Method for Sw\n");
02621 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
02622 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
02623 getPcFunction(&fPc,&dfPc,&fPcInv);
02624 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
02625 getPBoundaryCondition(m_fDP,m_fFluxBC);
02626 Function3D *fTransBC = getTransportForSystemBC(getFluxFunction()->numComponents());
02627 VecDouble &vK = getPermeabilityAtCells();
02628 VecDouble &vPor = getPorosityAtCells();
02629 LinearSolver *solver =getLinearSolver();
02630 m_pDiff = new CompDiffusiveStep(getMesh(),*m_fDP,*fTransBC,vK,vPor,dynamic_cast<Function1D&>(*fMobGrav),dynamic_cast<Function1D&>(*dfPc),*getFlashCompositional(),*solver);
02631
02632 m_pDiff->setTransport(*getTransportModule());
02633 m_pDiff->setDynamic(*getDynamicModule());
02634 break;
02635 }
02636 case 2:
02637 {
02638 print(section,"Double Porosity Diffusive Step\n");
02639 print(section,"\t(Celestin Stuff)\n");
02640
02641 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
02642 getPcFunction(&fPc,&dfPc,&fPcInv);
02643
02644 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
02645 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
02646
02647 m_pDiff=new DoublePorosityDiff(getMesh(),
02648 getPermeabilityAtCells(),
02649 getPorosityAtCells(),
02650 dynamic_cast<Function1D&>(*fMobGrav),
02651 dynamic_cast<Function1D&>(*dfPc),
02652 dynamic_cast<Function1D&>(*fPc),
02653 getBlockMatrix());
02654 break;
02655 }
02656 case 3:
02657 {
02658 print(section,"Diffusive Step for Biphasic Flow (see Douglas at.al)\n");
02659 print(section,"Note: No flux enter or leave the reservoir due to the\n\
02660 capillary pressure gradient\n");
02661
02662 GeneralFunctionInterface *fPc,*dfPc,*fPcInv;
02663 getPcFunction(&fPc,&dfPc,&fPcInv);
02664
02665 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav;
02666 getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav);
02667 m_pDiff=new BiphasicDiff(getMesh(),
02668 getPermeabilityAtCells(),
02669 getPorosityAtCells(),
02670 dynamic_cast<Function1D&>(*fMobGrav),
02671 dynamic_cast<Function1D&>(*dfPc),
02672 *getLinearSolver());
02673
02674 break;
02675 }
02676
02677 default:
02678 throw new Exception("Invalid Option for DIFFUSIVE_STEP");
02679 }
02680 return m_pDiff;
02681 }
02692 void MainApp::getMobilities(GeneralFunctionInterface* &fMob,GeneralFunctionInterface* &fMobT,GeneralFunctionInterface* &fFracFlux,GeneralFunctionInterface* &DfFracFlux,GeneralFunctionInterface* &fMobGrav,GeneralFunctionInterface* &DfMobGrav)
02693 {
02694 std::string section="Mobilities";
02695 if (!m_fFracFlux)
02696 {
02697 print(section,"Relative Mobility:\n\t");
02698
02699 switch(configFile.getInt("RELATIVE_MOBILITY_FUNCTION"))
02700 {
02701 case 1:
02702 {
02703 double v1 = configFile.getDouble("S1_VISCOSITY");
02704 double v2 = configFile.getDouble("S2_VISCOSITY");
02705 double Sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
02706 double Sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
02707 m_fMob= new FLinear(1.0/v1,0);
02708 m_fMobT=new FLinear((1.0/v1 - 1.0/v2),1.0/v2);
02709 m_fFracFlux = new FFractionalLinearMobility(v1,v2);
02710 m_DfFracFlux = new FDFractionalLinearMobility(v1,v2);
02711 m_DfMobGrav= new DFLinearMobilityProduct(v1,v2,Sr1,Sr2);
02712 m_fMobGrav= new FLinearMobilityProduct(v1,v2,Sr1,Sr2);
02713 break;
02714
02715 }
02716 case 2:
02717 {
02718 double v1 = configFile.getDouble("S1_VISCOSITY");
02719 double v2 = configFile.getDouble("S2_VISCOSITY");
02720 double Sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
02721 double Sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
02722 m_fMob = NULL;
02723 m_fFracFlux=new FChaventMobility(v1,v2,Sr1,Sr2);
02724 m_DfFracFlux=new DFChaventMobility(v1,v2,Sr1,Sr2);
02725 m_fMobT=new FBucleyLeverettTotalMob(v1,v2,Sr1,Sr2);
02726 m_fMobGrav=new FBucleyLeverettGravityMob(v1,v2,Sr1,Sr2);
02727 m_DfMobGrav=new DFBucleyLeverettGravityMob(v1,v2,Sr1,Sr2);
02728
02729
02730 print(section,"\
02731 mob1(S1)=(S1-Sr1)^2/(1-Sr1)^2\n\t\
02732 mob2(S1)=(1-S1/(1-Sr2))^2;\n\t\
02733 Residual Saturation (Sr1) = %g\n\t\
02734 Residual Saturation (Sr2) = %g\n\t\
02735 Viscosity of Phase 1 (v1) = %g\n\t\
02736 Viscosity of Phase 2 (v2) = %g\n\t\
02737 Maximum Char Velocity %g\n\t",Sr1,Sr2,v1,v2,
02738 dynamic_cast<Function1D&>(*m_DfFracFlux).getMaxNorm(0,1));
02739 break;
02740
02741
02742 }
02743 case 3:
02744 {
02745 print(section,"mob1(S1) = mob2(S2) = 1\n");
02746 m_fMob= new FQuadratic(0,0,1);
02747 break;
02748 }
02749 case 4:
02750 {
02751 double v1 = configFile.getDouble("S1_VISCOSITY");
02752 double v2 = configFile.getDouble("S2_VISCOSITY");
02753 double Sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
02754 double Sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
02755 print(section,"Mobilities based on the relative permeabilities\n\
02756 obtained from Piri-Morteza experiments\n\t\
02757 Residual Saturation (Sr1) = %g\n\t\
02758 Residual Saturation (Sr2) = %g\n\t\
02759 Viscosity of Water (v1) = %g\n\t\
02760 Viscosity of CO2 (v2) = %g\n\t",Sr1,Sr2,v1,v2);
02761 Function1D *krw= new FKrwPiri();
02762 Function1D *krc= new FKrcPiri();
02763 Function1D *dkrc=new DFKrcPiri();
02764 Function1D *dkrw=new DFKrwPiri();
02765 m_fKr1= krw;
02766 m_fKr2= krc;
02767 m_fDKr2=dkrc;
02768 m_fDKr1=dkrw;
02769
02770
02771 m_fMob = NULL;
02772 m_fFracFlux =new FFracFluxFromRelK ( *krw,*krc,Sr1,Sr2,v1,v2);
02773 m_DfFracFlux=new DFFracFluxFromRelK(*krw,*krc,*dkrw,*dkrc,Sr1,Sr2,v1,v2);
02774 m_fMobT =new FMobTFromRelK(*krw,*krc,Sr1,Sr2,v1,v2);
02775 m_fMobGrav =new FFracGravFromRelK(*krw,*krc,Sr1,Sr2,v1,v2);
02776 m_DfMobGrav =new DFFracGravFromRelK(*krw,*krc,*dkrw,*dkrc,Sr1,Sr2,v1,v2);
02777
02778 break;
02779 }
02780 default:
02781 {
02782 throw new Exception("Relative mobilities are not yet defined for the case\n\tFLUX_FUNCTION=%d\n",configFile.getInt("FLUX_FUNCTION"));
02783 }
02784 }
02785 }
02786 fMob=m_fMob;
02787 fFracFlux=m_fFracFlux;
02788 DfFracFlux=m_DfFracFlux;
02789 fMobT=m_fMobT;
02790 fMobGrav=m_fMobGrav;
02791 DfMobGrav=m_DfMobGrav;
02792
02793
02794 }
02809 void MainApp::getKr(GeneralFunctionInterface* &fKr,GeneralFunctionInterface* &DfKr)
02810 {
02811 if (!m_fKr)
02812 {
02813 switch(configFile.getInt("RELATIVE_MOBILITY_FUNCTION"))
02814 {
02815
02816 case 4:
02817 {
02818 Function1D *krw= new FKrwPiri();
02819 Function1D *krc= new FKrcPiri();
02820 Function1D *dkrc=new DFKrcPiri();
02821 Function1D *dkrw=new DFKrwPiri();
02822 m_fKr = new VectorFunction(krw,krc,true);
02823 m_fDKr=new VectorFunction(dkrw,dkrc,true);
02824 break;
02825 }
02826 default:
02827 m_fKr=NULL;
02828 m_fDKr=NULL;
02829
02830 }
02831 }
02832 fKr=m_fKr;
02833 DfKr=m_fDKr;
02834
02835 return;
02836 }
02837
02838
02839
02840
02851 void MainApp::mb_getMobilities(GeneralFunctionInterface* &fMob,GeneralFunctionInterface* &fMobT,GeneralFunctionInterface* &fFracFlux,GeneralFunctionInterface* &DfFracFlux,GeneralFunctionInterface* &fMobGrav,GeneralFunctionInterface* &DfMobGrav)
02852 {
02853 std::string section="Mobilities";
02854 if (!m_mb_fMob)
02855 {
02856 print(section,"Relative Mobility:\n\t");
02857
02858 switch(configFile.getInt("MB_RELATIVE_MOBILITY_FUNCTION"))
02859 {
02860 case 1:
02861 {
02862 double v1 = configFile.getDouble("S1_VISCOSITY");
02863 double v2 = configFile.getDouble("S2_VISCOSITY");
02864 m_mb_fMob= new FLinear(1.0/v1,0);
02865 m_mb_fMobT=new FLinear((1.0/v1 - 1.0/v2),1.0/v2);
02866 m_mb_fFracFlux = new FFractionalLinearMobility(v1,v2);
02867 m_mb_DfFracFlux = new FDFractionalLinearMobility(v1,v2);
02868 m_DfMobGrav=NULL;
02869 m_mb_fMobGrav=NULL;
02870 break;
02871
02872 }
02873 case 2:
02874 {
02875 double v1 = configFile.getDouble("S1_VISCOSITY");
02876 double v2 = configFile.getDouble("S2_VISCOSITY");
02877 double Sr1= configFile.getDouble("MB_PHASE1_RESIDUAL_SATURATION",0);
02878 double Sr2= configFile.getDouble("MB_PHASE2_RESIDUAL_SATURATION",0);
02879 print(section,"\
02880 mob1(S1)=(S1-Sr1)^2/(1-Sr1)^2\n\t\
02881 mob2(S1)=(1-S1/(1-Sr2))^2;\n\t\
02882 Residual Saturation (Sr1) = %g\n\t\
02883 Residual Saturation (Sr2) = %g\n\t\
02884 Viscosity of Phase 1 (v1) = %g\n\t \
02885 Viscosity of Phase 2 (v2) = %g\n",Sr1,Sr2,v1,v2);
02886 m_mb_fMob = NULL;
02887 m_mb_fFracFlux=new FChaventMobility(v1,v2,Sr1,Sr2);
02888 m_mb_DfFracFlux=new DFChaventMobility(v1,v2,Sr1,Sr2);
02889
02890 m_mb_fMobT=new FBucleyLeverettTotalMob(v1,v2,Sr1,Sr2);
02891 m_mb_fMobGrav=new FBucleyLeverettGravityMob(v1,v2,Sr1,Sr2);
02892 m_mb_DfMobGrav=new DFBucleyLeverettGravityMob(v1,v2,Sr1,Sr2);
02893 break;
02894
02895
02896 }
02897 case 3:
02898 {
02899 double v1 = configFile.getDouble("S1_VISCOSITY");
02900 double v2 = configFile.getDouble("S2_VISCOSITY");
02901 double Sr1= configFile.getDouble("MB_PHASE1_RESIDUAL_SATURATION");
02902 double Sr2= configFile.getDouble("MB_PHASE2_RESIDUAL_SATURATION");
02903 m_mb_fMob = new FSquareMob (v1, v2, Sr1);
02904 m_mb_fMobT = new FSquareMobT (v1, v2, Sr1, Sr2);
02905 m_mb_fFracFlux=NULL;
02906 m_mb_DfFracFlux=NULL;
02907 m_mb_fMobGrav=NULL;
02908 m_mb_DfMobGrav=NULL;
02909 break;
02910
02911 }
02912 case 4:
02913 {
02914 print(section,"mob1(S1) = mob2(S2) = 1\n");
02915 m_mb_fMob= new FQuadratic(0,0,1);
02916 break;
02917 }
02918 default:
02919 {
02920 throw new Exception("Relative mobilities are not yet defined for the case\n\tFLUX_FUNCTION=%d\n",configFile.getInt("FLUX_FUNCTION"));
02921 }
02922 }
02923 }
02924 fMob=m_mb_fMob;
02925 fFracFlux=m_mb_fFracFlux;
02926 DfFracFlux=m_mb_DfFracFlux;
02927 fMobT=m_mb_fMobT;
02928 fMobGrav=m_mb_fMobGrav;
02929 DfMobGrav=m_mb_DfMobGrav;
02930 }
02931
02932
02933
02934
02935 void MainApp::getPcFunction(GeneralFunctionInterface **fpc,GeneralFunctionInterface **dfpc,GeneralFunctionInterface **fpcInv )
02936 {
02937 std::string section="Pc";
02938 print(section,"Capillary Pressure:");
02939 if (!m_fPc)
02940 {
02941 switch(configFile.getInt("CAPILLARY_PRESSURE"))
02942 {
02943 case 0:
02944 {
02945 if (configFile.getInt("DIFFUSIVE_STEP",0)==0)
02946 print(section,"\n\tNull Cappillar Pressure\n");
02947 else {
02948 throw new Exception("WARNING: Null cappillar pressure will result in a non-singular system for the diffusive step\n");
02949 }
02950 m_fPc = new FPcLinear(0,0);
02951 m_dfPc = new FPcLinear(0,0);
02952 m_fPcInv=new FPcLinear(0,0);
02953 break;
02954 }
02955 case 1:
02956 {
02957 double d1 = configFile.getDouble("CAPILLARY_PRESSURE_ARG_1");
02958 double d2 = configFile.getDouble("CAPILLARY_PRESSURE_ARG_2");
02959 m_fPc = new FPcLinear(d1,d2);
02960 m_dfPc = new FLinear(0,d2-d1);
02961 m_fPcInv = new FLinear(1.0/(d2-d1),-d1/(d2-d1));
02962 print(section,"\tCapillar Pressure: Linear from %g to %g\n",d1,d2);
02963 break;
02964 }
02965 case 2:
02966 {
02967 double Sr1= configFile.getDouble("PHASE1_RESIDUAL_SATURATION",0);
02968 double Sr2= configFile.getDouble("PHASE2_RESIDUAL_SATURATION",0);
02969 double d1 = configFile.getDouble("CAPILLARY_PRESSURE_ARG_1");
02970 double d2 = configFile.getDouble("CAPILLARY_PRESSURE_ARG_2",INFINITY);
02971 print(section,"\
02972 \n\tCapillar Pressure: Pc(s) = eta ((s-Sr1)^-2 - tau (1-s)^-2)\
02973 \n\twhere\
02974 \n\teta=%g\
02975 \n\ttau=Sr2^2 * (1-Sr2-Sr1)^2\
02976 \n\tSr1=%g\
02977 \n\tSr2=%g\n\
02978 \n\tPressure Threeshold:%g\n",d1,Sr1,Sr2,d2);
02979 m_fPc = new FPcSquare(d1,Sr1,Sr2,d2);
02980 m_dfPc = new DFPcSquare(d1,Sr1,Sr2,d2);
02981 m_fPcInv = NULL;
02982 break;
02983 }
02984 case 3:
02985 {
02986 double gamma = configFile.getDouble("CAPILLARY_PRESSURE_ARG_1");
02987 double theta = configFile.getDouble("CAPILLARY_PRESSURE_ARG_2");
02988 m_fPc = new FPcFracture(gamma,theta);
02989 m_dfPc = new DFPcFracture(gamma,theta);
02990 m_fPcInv = new FPcInvFracture(gamma,theta);
02991 print(section,"\tCappillar Pressure: Fracture from %g to %g\n",gamma,theta );
02992 break;
02993 }
02994
02995 default:
02996 throw new Exception("Invalid option for CAPILLARY_PRESSURE");
02997
02998 }
02999 }
03000 *fpc = m_fPc;
03001 *dfpc = m_dfPc;
03002 *fpcInv = m_fPcInv;
03003 }
03004
03005 void MainApp::mb_getPcFunction(GeneralFunctionInterface **mb_fpc, GeneralFunctionInterface **mb_dfpc, GeneralFunctionInterface **mb_fpcInv)
03006 {
03007 string section = "MB_PC";
03008 if (!m_mb_fPc)
03009 {
03010 switch (configFile.getInt("MB_CAPILLARY_PRESSURE"))
03011 {
03012 case 0:
03013 {
03014 if (configFile.getInt("DIFFUSIVE_STEP",0)==0)
03015 print(section,"\n\tNull Capillar Pressure\n");
03016 else {
03017 throw new Exception("WARNING: Null capillar pressure will result in a non-singular system for the diffusive step\n");
03018 }
03019 m_mb_fPc = new FPcLinear(0,0);
03020 m_mb_dfPc = new FPcLinear(0,0);
03021 m_mb_fPcInv=new FPcLinear(0,0);
03022 break;
03023 }
03024 case 1:
03025 {
03026 double d1 = configFile.getDouble("MB_CAPILLARY_PRESSURE_ARG_1");
03027 double d2 = configFile.getDouble("MB_CAPILLARY_PRESSURE_ARG_2");
03028 m_mb_fPc = new FPcLinear(d1,d2);
03029 m_mb_dfPc = new FLinear(0,d2-d1);
03030 m_mb_fPcInv = new FLinear(1.0/(d2-d1),-d1/(d2-d1));
03031 print(section,"\tCapillar Pressure: Linear from %g to %g\n",d1,d2);
03032 break;
03033 }
03034 case 2:
03035 {
03036 double Sr1 = configFile.getDouble("MB_PHASE1_RESIDUAL_SATURATION",0);
03037 double Sr2 = configFile.getDouble("MB_PHASE2_RESIDUAL_SATURATION",0);
03038 double maxS1 = 1.0-Sr2;
03039 double d1 = configFile.getDouble("MB_CAPILLARY_PRESSURE_ARG_1");
03040
03041 double S_max = configFile.getDouble("MB_CAPILLARY_PRESSURE_ARG_2");
03042 unsigned N_subints = configFile.getUnsigned("MB_CAPILLARY_PRESSURE_ARG_3");
03043 double BS_tol=configFile.getDouble("BS_METHOD_TOLERANCE",1.0e-6);
03044 double P_max_Flash = configFile.getDouble("MB_CAPILLARY_PRESSURE_ARG_4", INFINITY);
03045
03046 print(section,"\
03047 \n\tCapillar Pressure: Pc(s) = eta ((s-Sr1)^-2 - tau (1-s)^-2)\
03048 \n\twhere\
03049 \n\teta=%g\
03050 \n\ttau=Sr2^2 * (1-Sr2-Sr1)^2\
03051 \n\tSr1=%g\
03052 \n\tSr2=%g\
03053 \n\tS_max=%g\
03054 \n\tN_subints=%d\
03055 \n\tBS_tol=%g\
03056 \n\tP_max_Flash=%g\
03057 \n",d1,Sr1,Sr2,S_max,N_subints,BS_tol,P_max_Flash);
03058
03059 m_mb_fPc = new FPcSquare(d1,Sr1,Sr2,P_max_Flash);
03060 m_mb_dfPc = new DFPcSquare(d1,Sr1,Sr2,P_max_Flash);
03061 m_mb_fPcInv = new InvCapPress(dynamic_cast<Function1D&>(*m_mb_fPc), Sr1, maxS1, S_max, N_subints, BS_tol);
03062 break;
03063 }
03064 default:
03065 throw new Exception ("Invalid option for Matrix Block Capillary Pressure");
03066 }
03067 }
03068
03069 *mb_fpc = m_mb_fPc;
03070 *mb_dfpc= m_mb_dfPc;
03071 *mb_fpcInv = m_mb_fPcInv;
03072
03073 }
03074
03075 void MainApp::printLog(std::ostream &out)
03076 {
03077 printBeginSection("Triangulation",out);
03078 print_log_entry("Mesh",out);
03079
03080 printBeginSection("Porous Media Properties",out);
03081 print_log_entry("PERMEABILITY",out);
03082 print_log_entry("POROSITY",out);
03083
03084 printBeginSection("Boundary Conditions",out);
03085 print_log_entry("TransportBC",out);
03086 print_log_entry("PBC",out);
03087
03088 printBeginSection("Initial Conditions",out);
03089 print_log_entry("PressureIC",out);
03090 print_log_entry("IC",out);
03091
03092
03093 if (! log["Dynamic"].empty())
03094 {
03095 printBeginSection("Dynamic Module",out);
03096 out << log["Dynamic"];
03097 out << log["LINEAR_SOLVER"];
03098 }
03099
03100 if (! log["Transport"].empty())
03101 {
03102 printBeginSection("Transport Method",out);
03103 out << log["Transport"];
03104
03105 }
03106
03107 if (! log["Flash"].empty())
03108 {
03109 printBeginSection("Flash",out);
03110 out << log["Flash"];
03111 }
03112
03113
03114 if (! log["Wells"].empty())
03115 {
03116 printBeginSection("Wells",out);
03117 out << log["Wells"];
03118 }
03119
03120
03121
03122 printBeginSection("Transport Convective Flux",out);
03123 print_log_entry("Flux",out);
03124
03125
03126 printBeginSection("Phase Dynamic",out);
03127 print_log_entry("Pc",out);
03128 print_log_entry("Mobilities",out);
03129
03130 printBeginSection("Sequence Algorithm",out);
03131 print_log_entry("Sequencer",out);
03132
03133 printBeginSection("End Log",out);
03134
03135 out << std::endl;
03136 }
03137
03138
03139 void MainApp::print_log_entry(std::string section,std::ostream &out)
03140 {
03141 if (!log[section].empty())
03142 out << log[section];
03143 }
03144
03148 MonitorWells& MainApp::getMonitorWells()
03149 {
03150 if (!m_pMonitorWells)
03151 {
03152 m_pMonitorWells = new MonitorWells();
03153
03154 std::string keyX="MONITOR_WELL_X_";
03155 std::string keyY="MONITOR_WELL_Y_";
03156 std::string keyZ="MONITOR_WELL_Z_";
03157
03158 unsigned i=0;
03159 std::string strI = StringFunctions::to_string(i++);
03160 MonitorWells::_NodeWell well;
03161 OrthoMesh &mesh = getMesh();
03162
03163 while (configFile.isDefined(keyX+strI))
03164 {
03165 well.P[0]=configFile.getDouble(keyX+strI);
03166 well.P[1]=configFile.getDouble(keyY+strI);
03167 well.P[2]=configFile.getDouble(keyZ+strI);
03168 well.cell_index = mesh.getCellAt(well.P)->index();
03169 m_pMonitorWells->getWells().push_back(well);
03170 }
03171 m_pMonitorWells->openFileName(getMonitorWellsOutputFileName());
03172 }
03173 return *m_pMonitorWells;
03174 }
03175
03176
03177
03178
03179 FlashCompositional* MainApp::getFlashCompositional()
03180 {
03181 FlashBase *flash = getFlash();
03182
03183 if (dynamic_cast<FlashCompositional*>(flash) != NULL)
03184 return dynamic_cast<FlashCompositional*>( flash);
03185 else
03186 throw new Exception("Flash Module Selected is expected to be a Compositional Flash");
03187 }
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
