#include <mainapp.h>
Definition at line 30 of file mainapp.h.
| MainApp::MainApp | ( | ) |
Definition at line 100 of file mainapp.cpp.
00101 { 00102 pOut = &(std::cout); 00103 //Initialize all the pointers for all the classes 00104 //as NULL 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 }
| MainApp::~MainApp | ( | ) |
Definition at line 122 of file mainapp.cpp.
00123 { 00124 //Desallocate all the pointers; 00125 //Remeber that the destructor can be called from a error exception. 00126 //In this cases not all the pointers have valid objects, some of them 00127 00128 //may have its initial value (NULL) (see the constructor of this class) 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 }
| void MainApp::abortEvent | ( | ) |
Definition at line 2559 of file mainapp.cpp.
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 }
| void MainApp::adjustInitialConditionForCompressibleModel | ( | DynamicBase & | dyn, | |
| ConservativeMethodForSystem & | trans, | |||
| FlashCompositional & | flash | |||
| ) | [protected] |
If the problem is compressible, the volume of the mixture changes with pressure. So we must adjust the initial condition of the transport such that the volume of all the mixture fills all the porous media.This is for consistence purpose. We do that by changing the total amount of moles trying to preserve the proportion of moles of each component.
| dyn | Dynamic module | |
| trans | Transport Module | |
| flash | Flash Module |
Definition at line 2360 of file mainapp.cpp.
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 try 02371 { 02372 MixedHybridCompositionalSimple *pDyn = dynamic_cast<MixedHybridCompositionalSimple*>(&dyn); 02373 pDyn->findInitialPressure(100,1.e-7); 02374 } 02375 catch (std::bad_cast* ) 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 //printf("%g %g %g\n",factor,phasesVol(0),phasesVol(1)); 02387 //flashData.print(); 02388 for (unsigned i=0;i<nComps;i++) 02389 trans.getSolutionAtCells()(cell,i)*=factor; 02390 } 02391 flash.execute(); 02392 //Get fixed pressure condition 02393 Function3D *fDP,*fFlux; 02394 getPBoundaryCondition(fDP,fFlux); 02395 trans.updateFixedPressureBC(fDP,flash); 02396 02397 02398 //Now set the 02399 02400 // for (unsigned cell=0;cell<nCells;cell++) 02401 // { 02402 // flash.flash(cell,flashData); 02403 // flash.getPhasesVolume(dyn.getPressureAtCells()(cell),flashData,phasesVol); 02404 // printf("%d = %g = %g\n",cell,(phasesVol(0)+phasesVol(1)),volCell); 02405 // } 02406 02407 02408 }
| void MainApp::appendSuffixToFileNames | ( | std::string | suffix | ) | [inline, protected] |
| void MainApp::coarseMesh | ( | unsigned | nX, | |
| unsigned | nY, | |||
| unsigned | nZ, | |||
| string | fileIn, | |||
| string | fileOut | |||
| ) |
Definition at line 2252 of file mainapp.cpp.
02253 { 02254 try{ 02255 //print arguments 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 //Initialize the HDF5 reader and writer. 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 }
| void MainApp::enablePrint | ( | bool | b | ) | [inline, protected] |
| void MainApp::execute | ( | std::string | fileName, | |
| bool | unitests = false | |||
| ) |
Execute the program
Definition at line 682 of file mainapp.cpp.
00683 { 00684 00685 try { 00686 //Read the config file. 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 //Initialize HDF5 writer with the triangulation. 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: //Case to solve the biphasic incompressible flow 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 /* Now if this is a compressible Model the initial condition of the transport 00774 must be updated such that all mixture could have the volume filled buy all 00775 the pores. We do that buy changing the number of total moles of all the mixture 00776 preserving the same proportion of moles in each component previously 00777 defined buy the initial conditions of the transport equations 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 //if (!hasSystemOfTransportEquations()) 00788 // throw new Exception("SEQUENCE_ALGHORITHM=3 was designed just for system of transport equations"); 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; //Just to avoid compiling warnings 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 }
| BlockMatrixModule & MainApp::getBlockMatrix | ( | ) | [protected] |
Definition at line 2571 of file mainapp.cpp.
02572 { 02573 //Get the mobilities function 02574 GeneralFunctionInterface *fMob,*fFracF,*fMobT,*fMobGrav,*DfFracF,*DfMobGrav; 02575 mb_getMobilities(fMob,fMobT,fFracF,DfFracF,fMobGrav,DfMobGrav); 02576 //Get the capillary pressure functions 02577 GeneralFunctionInterface *fPc,*dfPc,*fPcInv; 02578 mb_getPcFunction(&fPc,&dfPc,&fPcInv); 02579 //pcIc = The initial pc from the fracture system 02580 VecDouble pcIC(getMesh().numCells()); 02581 // { 02582 // GeneralFunctionInterface *fPc,*dfPc,*fPcInv; 02583 // getPcFunction(&fPc,&dfPc,&fPcInv); 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 }
| ConstVelocityModule * MainApp::getConstVelocityDynamicModule | ( | ) | [protected] |
Definition at line 1903 of file mainapp.cpp.
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; //Just to avoid warnings 01944 } 01945 }
| DiffusiveStep * MainApp::getDiffusiveStep | ( | ) | [protected] |
Definition at line 2599 of file mainapp.cpp.
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 //Get Pc function 02641 GeneralFunctionInterface *fPc,*dfPc,*fPcInv; 02642 getPcFunction(&fPc,&dfPc,&fPcInv); 02643 //Get the mobilities 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 //Get Pc function 02662 GeneralFunctionInterface *fPc,*dfPc,*fPcInv; 02663 getPcFunction(&fPc,&dfPc,&fPcInv); 02664 //Get the mobilities 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 }
| Point3D MainApp::getDomainDimensions | ( | ) | [protected] |
Definition at line 1833 of file mainapp.cpp.
01833 { 01834 return Point3D(configFile.getDouble("DIM_X"), 01835 configFile.getDouble("DIM_Y"), 01836 configFile.getDouble("DIM_Z")); 01837 }
| Point3D MainApp::getDX | ( | ) | [protected] |
Get the discretizations steps in eahc direction and store them in the parameter p
| p1 | (Output) Point in 3D containing the discretization steps in each direction |
Definition at line 2000 of file mainapp.cpp.
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 }
| DynamicBase * MainApp::getDynamicModule | ( | ) | [protected] |
Load the dynamic module
Create the dynamic module
| tria |
Definition at line 872 of file mainapp.cpp.
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 //Permeability (Todo: Functions to get permeability) 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 //Permeability (Todo: Functions to get permeability) 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 //Permeability (Todo: Functions to get permeability) 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 //m_fK = getK(); 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 //Permeability (Todo: Functions to get permeability) 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 //Permeability (Todo: Functions to get permeability) 01021 m_pDynamicModule= new PoissonFEM(getMesh(),getPermeabilityAtCells(), 01022 *m_fDP,*m_fFluxBC, 01023 *getLinearSolver(), 01024 configFile.getUnsigned("DYNAMIC_DEBUG_LEVEL")); 01025 //,*m_fGravMobSrc, 01026 // configFile.getDouble("TOLERANCE"), 01027 // configFile.getInt("NUM_ITERATIONS"), 01028 // configFile.getUnsigned("DYNAMIC_DEBUG_LEVEL")); 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 //Permeability (Todo: Functions to get permeability) 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 //Permeability (Todo: Functions to get permeability) 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 }
| double MainApp::getDynamicTimeStep | ( | ) | [protected] |
Definition at line 2303 of file mainapp.cpp.
02304 { 02305 double nIt = configFile.getDouble("N_DYN_ITERATIONS"); 02306 double time = configFile.getDouble("END_TIME"); 02307 return time/nIt; 02308 }
| FlashBase * MainApp::getFlash | ( | ) | [protected] |
Definition at line 508 of file mainapp.cpp.
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 //double v2 =configFile.getDouble("S1_VISCOSITY"); 00614 //double v1 =configFile.getDouble("S1_VISCOSITY"); 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 //double v2 =configFile.getDouble("S1_VISCOSITY"); 00666 //double v1 =configFile.getDouble("S1_VISCOSITY"); 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 }
| FlashCompositional * MainApp::getFlashCompositional | ( | ) | [protected] |
Definition at line 3179 of file mainapp.cpp.
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 }
| FaceFluxFunction * MainApp::getFluxFunction | ( | ) | [protected] |
Get the faceflux object specifing the flux term of the conservative equations.
Definition at line 1489 of file mainapp.cpp.
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 //return new FaceFluxWithGravity(getMesh().numFaces(),*new FLinear(0,0),*pGravMob,*m_pCK,-(pw-po)*g); 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 //m_flux= new FaceFluxWithIndependentEquations(getMesh(),new FChaventMobility(1,20),new FLinear(1,0),true); 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 }
| std::string MainApp::getHDF5FileName | ( | ) | [protected] |
Definition at line 2523 of file mainapp.cpp.
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 }
| void MainApp::getKr | ( | GeneralFunctionInterface *& | fKr, | |
| GeneralFunctionInterface *& | DfKr | |||
| ) | [protected] |
Return the relative permeabilities. This function has a strong connection with the MainApp::getMobilities in the sense that they read the same keys from the config file and the mobilities are defined based on the values of the relative permeabilities; Note that depending of the option of the key RELATIVE_MOBILITY_FUNCTION the relative permeability can be not defined and in this case the output arguments fKr and DfKr return NULL. The relative mobilities returned are vector functions where the first component is the relative permeability of phase 1, second component for the relative permeability of phase 2 and so on fKr To containt the relative permeabilities DfKr To contain the derivative of the relative permeabilities
Definition at line 2809 of file mainapp.cpp.
02810 { 02811 if (!m_fKr) 02812 { 02813 switch(configFile.getInt("RELATIVE_MOBILITY_FUNCTION"/*,2)*/)) 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 }
| LinearSolver * MainApp::getLinearSolver | ( | ) | [protected] |
Definition at line 2311 of file mainapp.cpp.
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 }
| OrthoMesh & MainApp::getMB_Mesh | ( | ) | [protected] |
Definition at line 409 of file mainapp.cpp.
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 }
| OrthoMesh & MainApp::getMesh | ( | ) | [protected] |
Definition at line 381 of file mainapp.cpp.
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 //delete all wells that are not supposed to be represented as holes in the domain 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 }
| void MainApp::getMobilities | ( | GeneralFunctionInterface *& | fMob, | |
| GeneralFunctionInterface *& | fMobT, | |||
| GeneralFunctionInterface *& | fFracFlux, | |||
| GeneralFunctionInterface *& | DfFracFlux, | |||
| GeneralFunctionInterface *& | fMobGrav, | |||
| GeneralFunctionInterface *& | DfMobGrav | |||
| ) | [protected] |
Get Mobilities
| fMob | To contain the phase mobility | |
| fMobT | Total phases mobility | |
| fFracFlux | Fractional mobility function | |
| DfFracFlux | Derivate of the fractional mobility function | |
| fMobGrav | The gravity mobility, defined as the product of the water, oil relative mobilities and the total mobility (i.e LambdaW LambdaT LambdaO | |
| DfMobGrav | The derivative of fMobGrav |
Definition at line 2692 of file mainapp.cpp.
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"/*,2)*/)) 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 }
| MonitorWells & MainApp::getMonitorWells | ( | ) | [protected] |
Get monitor wells
Definition at line 3148 of file mainapp.cpp.
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 }
| std::string MainApp::getMonitorWellsOutputFileName | ( | ) | [protected] |
Definition at line 2547 of file mainapp.cpp.
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 }
| std::vector< unsigned > MainApp::getNElements | ( | ) | [protected] |
Definition at line 2140 of file mainapp.cpp.
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 }
| Point3D MainApp::getP0 | ( | ) | [protected] |
Definition at line 1981 of file mainapp.cpp.
01982 { 01983 return Point3D(0.0,0.0,0.0); 01984 }
| Point3D MainApp::getP1 | ( | ) | [protected] |
Definition at line 1989 of file mainapp.cpp.
01990 { 01991 return Point3D(configFile.getDouble("DIM_X"), 01992 configFile.getDouble("DIM_Y"), 01993 configFile.getDouble("DIM_Z")); 01994 }
| void MainApp::getPBoundaryCondition | ( | Function3D *& | fDP, | |
| Function3D *& | fFlux | |||
| ) | [protected] |
GetBoundary Conditions for the pressure using the two output parameters.
| fDP | It contains the dirichlet conditions for the pressure. | |
| fFlux | It contains the conditions for the prescribed flux. Actually such function define the normal component of the velocity along the boundary |
Definition at line 1124 of file mainapp.cpp.
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"), /* Now if this is a compressible Model the initial condition of the transport 01142 must be updated such that all mixture could have the volume filled buy all 01143 the pores. We do that buy changing the number of total moles of all the mixture 01144 preserving the same proportion of moles in each component previously 01145 defined buy the initial conditions of the transport equations 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);//Well Inj 01281 Function3D *fAux1 = new FCylinderRegion(C0,radius,height,vI);//Well Inj 01282 C0[0]=getMesh().getQ()[0]; 01283 C0[1]=getMesh().getQ()[1]; 01284 fDP = new FCylinderRegion(C0,radius,height,pO);//Well Prod 01285 Function3D *fAux2 = new FCylinderRegion(C0,radius,height,pO);//Well Prod 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 }
| void MainApp::getPcFunction | ( | GeneralFunctionInterface ** | fpc, | |
| GeneralFunctionInterface ** | dfpc, | |||
| GeneralFunctionInterface ** | fpcInv | |||
| ) | [protected] |
Definition at line 2935 of file mainapp.cpp.
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;//(1,2,3,4,5,6,7,8,9);//PcInvFunc(d1,s1,s2,s3,ipcm1,ipcm2,ipcm3,Sr1,Sr2); // J. Douglas et. al m_fPc = 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 }
| VecDouble & MainApp::getPermeabilityAtCells | ( | ) | [protected] |
Definition at line 1866 of file mainapp.cpp.
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 //plot this field 01877 VecDouble vV(getMesh().numVertices()); 01878 getMesh().projectCentralValuesAtVertices(m_cK,vV); 01879 HDF5OrthoWriter::getHDF5OrthoWriter().writeScalarField(vV,"K"); 01880 } 01881 return m_cK; 01882 01883 }
| double MainApp::getPoisson | ( | ) | [protected] |
Get Poisson elastic modulus
Definition at line 2294 of file mainapp.cpp.
02295 { 02296 double result = configFile.getDouble("POISSON_COEFFICIENT"); 02297 print("POISSON","POISSON COEFFICIENT = %g\n",result); 02298 return result; 02299 }
| VecDouble & MainApp::getPorosityAtCells | ( | ) | [protected] |
Definition at line 1842 of file mainapp.cpp.
| Function3D * MainApp::getPorousMediaProperty | ( | std::string | key | ) | [protected] |
Definition at line 1410 of file mainapp.cpp.
01411 { 01412 std::string section=key; 01413 //Get the names of the keys used as arguments 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 }
| FixedValueCondition & MainApp::getPressureFixedConditions | ( | ) | [protected] |
Definition at line 2495 of file mainapp.cpp.
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++) //for all cells 02504 { 02505 //get barycenter 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 }
| Function3D * MainApp::getPressureIC | ( | ) | [protected] |
Get Initial condition for pressure
Definition at line 334 of file mainapp.cpp.
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; //The return value is here just to avoid compiling warnings 00354 break; 00355 } 00356 return NULL; 00357 }
| Sequencer * MainApp::getSequencer | ( | ) | [protected] |
Load the sequencer
Definition at line 359 of file mainapp.cpp.
00360 { 00361 00362 return new Sequencer(getMonitorWells()); 00363 }
| VecWellInfo MainApp::getSourceWells | ( | ) | [protected] |
Definition at line 1816 of file mainapp.cpp.
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 }
| Point3D MainApp::getSpatialDiscretizationSteps | ( | ) | [protected] |
Get the discretizations steps in eahc direction and store them in the parameter p
| p1 | (Output) Point in 3D containing the discretization steps in each direction |
Definition at line 1387 of file mainapp.cpp.
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 }
| Function3D * MainApp::getTransportBC | ( | std::string | compName | ) | [protected] |
Definition at line 150 of file mainapp.cpp.
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 // radius+=std::max(getMesh().getDX()[0],getMesh().getDX()[1]); 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 //Print a error message and quits the application 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 }
| FixedValueCondition & MainApp::getTransportFixedConditions | ( | ) | [protected] |
Definition at line 2445 of file mainapp.cpp.
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 //Now append fixed conditions for the wells 02482 //Get only the wells that are implemented as source terms 02483 VecWellInfo wells = getSourceWells(); 02484 02485 02486 02487 02488 02489 m_fixedC. addTransportFixedCondition(mesh,wells); 02490 firstTime=false; 02491 } 02492 return m_fixedC; 02493 }
| Function3D * MainApp::getTransportForSystemBC | ( | unsigned | nCmps | ) | [protected] |
Get boundary conditions for the system of transport equations
| flux |
Definition at line 2416 of file mainapp.cpp.
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 }
| Function3D * MainApp::getTransportForSystemIC | ( | unsigned | nCmps | ) | [protected] |
Get initial conditions for the system of transport equations
| flux |
Definition at line 2432 of file mainapp.cpp.
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 }
| Function3D * MainApp::getTransportIC | ( | std::string | compName | ) | [protected] |
Definition at line 260 of file mainapp.cpp.
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 //Print message error and quit the application 00324 throw new Exception("Invalid option for %s_INITIAL_VALUE",compName.c_str()); 00325 return NULL; //The return value is here just to avoid compiling warnings 00326 break; 00327 } 00328 return NULL; 00329 }
| TransportBase * MainApp::getTransportModule | ( | ) | [protected] |
Load the tranport module for scalar equations
Definition at line 427 of file mainapp.cpp.
00428 { 00429 std::string section="Transport"; 00430 if (m_pTransportModule) 00431 return m_pTransportModule; 00432 00433 00434 //Get Boundary conditions and initial conditions 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 //Print error message and quit the application 00500 throw new Exception("Invalid option for TRANSPORT_MODULE selection"); 00501 return NULL; //The return value is here just to avoid compiling warnings 00502 00503 } 00504 return m_pTransportModule; 00505 }
| void MainApp::getUBoundaryCondition | ( | Function3D *& | fDU, | |
| Function3D *& | fTensionBC | |||
| ) | [protected] |
GetBoundary Conditions for the displacement using the two output parameters.
| fDU | Contains the dirichlet conditions for the displacement. | |
| fTensionBC | Will contain the boundary conditions for the tensions. It is a vector function, each component represents the integral surface intS ( Ti * n) where Ti is the ith line of the tensor of tensions T and n represents the normal at the boundary |
Definition at line 1338 of file mainapp.cpp.
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 }
| VecWellInfo MainApp::getWells | ( | ) | [protected] |
Definition at line 1666 of file mainapp.cpp.
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 //Now adjust the boundaries of each well. 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 //Check the wells dont intersect the boundary of the domain 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 }
| void MainApp::getWellsPressureBC | ( | Function3D *& | fP, | |
| Function3D *& | fFlux | |||
| ) | [protected] |
Definition at line 2011 of file mainapp.cpp.
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 //Well with pressure condition 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 }
| Function3D & MainApp::getWellsSourceTerm | ( | ) | [protected] |
Definition at line 2060 of file mainapp.cpp.
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 }
| Function3D * MainApp::getWellTransportBC | ( | std::string | varName | ) | [protected] |
Definition at line 2091 of file mainapp.cpp.
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 }
| double MainApp::getYoung | ( | ) | [protected] |
Get Young elastic modulus
Definition at line 2287 of file mainapp.cpp.
02288 { 02289 double result = configFile.getDouble("YOUNG_COEFFICIENT"); 02290 print("YOUNG","YOUNG COEFFICIENT = %g\n",result); 02291 return result; 02292 }
| unsigned MainApp::hasSystemOfTransportEquations | ( | ) | [protected] |
Definition at line 1954 of file mainapp.cpp.
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 }
| bool MainApp::isCompositionalFlash | ( | ) | [protected] |
| void MainApp::mb_getMobilities | ( | GeneralFunctionInterface *& | fMob, | |
| GeneralFunctionInterface *& | fMobT, | |||
| GeneralFunctionInterface *& | fFracFlux, | |||
| GeneralFunctionInterface *& | DfFracFlux, | |||
| GeneralFunctionInterface *& | fMobGrav, | |||
| GeneralFunctionInterface *& | DfMobGrav | |||
| ) | [protected] |
Get Mobilities for matrix block: dynamic module fMob To contain the phase mobility fMobT Total phases mobility fFracFlux Fractional mobility function DfFracFlux Derivate of the fractional mobility function fMobGrav The gravity mobility, defined as the product of the water, oil relative mobilities and the total mobility (i.e LambdaW LambdaT LambdaO DfMobGrav The derivative of fMobGrav
Definition at line 2851 of file mainapp.cpp.
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"/*,2)*/)) 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; //Todo: Not implemented yet 02869 m_mb_fMobGrav=NULL; //Todo: Not implemented yet 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"/*,0*/); 02902 double Sr2= configFile.getDouble("MB_PHASE2_RESIDUAL_SATURATION"/*,0*/); 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 }
| void MainApp::mb_getPcFunction | ( | GeneralFunctionInterface ** | mb_fpc, | |
| GeneralFunctionInterface ** | mb_dfpc, | |||
| GeneralFunctionInterface ** | mb_fpcInv | |||
| ) | [protected] |
Definition at line 3005 of file mainapp.cpp.
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); // J. Douglas et. al 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 }
| VecDouble & MainApp::mb_getPermeabilityAtCells | ( | ) | [protected] |
Definition at line 1885 of file mainapp.cpp.
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 }
| VecDouble & MainApp::mb_getPorosityAtCells | ( | ) | [protected] |
Definition at line 1853 of file mainapp.cpp.
| void MainApp::print | ( | std::string | section, | |
| const char * | format, | |||
| ... | ||||
| ) | [protected] |
Definition at line 2121 of file mainapp.cpp.
| void MainApp::print_log_entry | ( | std::string | section, | |
| std::ostream & | out | |||
| ) | [protected] |
Definition at line 3139 of file mainapp.cpp.
| void MainApp::printBeginSection | ( | std::string | str, | |
| std::ostream & | out = std::cout | |||
| ) | [protected] |
Definition at line 1083 of file mainapp.cpp.
| void MainApp::printLog | ( | std::ostream & | out | ) | [protected] |
Definition at line 3075 of file mainapp.cpp.
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 }
| void MainApp::printMeshInfo | ( | ) | [protected] |
Print information about triangulation and other structures
Definition at line 1400 of file mainapp.cpp.
01401 { 01402 if (configFile.getInt("DEBUG_PRINT_TRIANGULATION")) 01403 getMesh().print(); 01404 01405 }
| void MainApp::printTriangulation | ( | ) | [protected] |
Definition at line 1097 of file mainapp.cpp.
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 }
| void MainApp::setOutput | ( | std::ostream & | out | ) | [inline, protected] |
| void MainApp::setTriangulation | ( | Triangulation< 3 > & | tria | ) | [protected] |
Definition at line 365 of file mainapp.cpp.
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 //OK generate the grid 00378 CO2GridGenerator::subdivided_hyper_rectangle_with_hole(tria,v,P1,P2,wells); 00379 }
| void MainApp::textcolor | ( | int | attr, | |
| int | fg, | |||
| int | bg | |||
| ) | [protected] |
Definition at line 1107 of file mainapp.cpp.
01108 { 01109 char command[13]; 01110 01111 /* Command is the control command to the terminal */ 01112 sprintf(command, "%c[%d;%d;%dm", 0x1B, attr, fg + 30, bg + 40); 01113 print("%s", command); 01114 }
ConfigFile MainApp::configFile [protected] |
Dictionary MainApp::log [protected] |
BlockMatrixModule* MainApp::m_blockMatrix [protected] |
bool MainApp::m_bPrint [protected] |
VecDouble MainApp::m_cK [protected] |
VecDouble MainApp::m_cPor [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_DfFracFlux [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_DfMobGrav [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_dfPc [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fDKr [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fDKr1 [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fDKr2 [protected] |
OPointer<Function3D> MainApp::m_fDP [protected] |
OPointer<Function3D> MainApp::m_fDU [protected] |
OPointer<Function3D> MainApp::m_fFluxBC [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fFracFlux [protected] |
Function3D* MainApp::m_fInit [protected] |
OPointer<Function3D> MainApp::m_fInitPressure [protected] |
FixedValueCondition MainApp::m_fixedC [protected] |
FixedValueCondition MainApp::m_fixedP [protected] |
Function3D* MainApp::m_fK [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fKr [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fKr1 [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fKr2 [protected] |
FaceFluxFunction* MainApp::m_flux [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fMob [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fMobGrav [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fMobT [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fPc [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_fPcInv [protected] |
Function3D* MainApp::m_fPor [protected] |
OPointer<Function3D> MainApp::m_fSource [protected] |
OPointer<Function3D> MainApp::m_fTensionBC [protected] |
Function3D* MainApp::m_fTransBC [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_DfFracFlux [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_DfMobGrav [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_dfPc [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fFracFlux [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fMob [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fMobGrav [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fMobT [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fPc [protected] |
OPointer<GeneralFunctionInterface> MainApp::m_mb_fPcInv [protected] |
VecDouble MainApp::m_mbcK [protected] |
VecDouble MainApp::m_mbcPor [protected] |
Function3D* MainApp::m_mbfK [protected] |
Function3D* MainApp::m_mbfPor [protected] |
OPointer<DiffusiveStep> MainApp::m_pDiff [protected] |
OPointer<DynamicBase> MainApp::m_pDynamicModule [protected] |
FlashBase* MainApp::m_pFlash [protected] |
OPointer<OrthoMesh> MainApp::m_pMB_Mesh [protected] |
OPointer<OrthoMesh> MainApp::m_pMesh [protected] |
OPointer<MonitorWells> MainApp::m_pMonitorWells [protected] |
OPointer<Function1D> MainApp::m_pRelMob [protected] |
Sequencer* MainApp::m_pSequencer [protected] |
OPointer<LinearSolver> MainApp::m_pSolver [protected] |
OPointer<TransportBase> MainApp::m_pTransportModule [protected] |
Function3D* MainApp::m_pVelFunction [protected] |
std::string MainApp::m_suffix [private] |
VecWellInfo MainApp::m_wells [protected] |
std::ostream* MainApp::pOut [private] |
1.6.3
