yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11614
[Branch ~yade-pkg/yade/git-trunk] Rev 3520: -move initialization() to 2PF
------------------------------------------------------------
revno: 3520
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Mon 2014-11-03 19:59:34 +0100
message:
-move initialization() to 2PF
modified:
pkg/pfv/TwoPhaseFlowEngine.cpp
pkg/pfv/TwoPhaseFlowEngine.hpp
pkg/pfv/UnsaturatedEngine.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/pfv/TwoPhaseFlowEngine.cpp'
--- pkg/pfv/TwoPhaseFlowEngine.cpp 2014-10-22 09:15:54 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp 2014-11-03 18:59:34 +0000
@@ -431,5 +431,77 @@
vtkfile.end_data();
}
+void TwoPhaseFlowEngine::initialization()
+{
+ scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene
+ setPositionsBuffer(true);//copy sphere positions in a buffer...
+ buildTriangulation(0.0,*solver);//create a triangulation and initialize pressure in the elements (connecting with W-reservoir), everything will be contained in "solver"
+// initializeCellIndex();//initialize cell index
+ computePoreThroatRadius();//save pore throat radius before drainage. Thomas, here you can also revert this to computePoreThroatCircleRadius().
+ computePoreBodyRadius();//save pore body radius before imbibition
+ computePoreBodyVolume();//save capillary volume of all cells, for fast calculating saturation
+ computeSolidLine();//save cell->info().solidLine[j][y]
+ initializeReservoirs();//initial pressure, reservoir flags and local pore saturation
+ solver->noCache = true;
+}
+
+void TwoPhaseFlowEngine::computeSolidLine()
+{
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = Tri.finite_cells_end();
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
+ for(int j=0; j<4; j++) {
+ solver->lineSolidPore(cell, j);
+ }
+ }
+ if(solver->debugOut) {cout<<"----computeSolidLine-----."<<endl;}
+}
+
+void TwoPhaseFlowEngine::initializeReservoirs()
+{
+ boundaryConditions(*solver);
+ solver->pressureChanged=true;
+ solver->reApplyBoundaryConditions();
+ ///keep boundingCells[2] as W-reservoir.
+ for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
+ (*it)->info().isWRes = true;
+ (*it)->info().isNWRes = false;
+ (*it)->info().saturation=1.0;
+ }
+ ///keep boundingCells[3] as NW-reservoir.
+ for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
+ (*it)->info().isNWRes = true;
+ (*it)->info().isWRes = false;
+ (*it)->info().saturation=0.0;
+ }
+
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ ///if we start from drainage
+ if(drainageFirst)
+ {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().Pcondition) continue;
+ cell->info().p()=bndCondValue[2];
+ cell->info().isWRes = true;
+ cell->info().isNWRes= false;
+ cell->info().saturation=1.0;
+ }
+ }
+ ///if we start from imbibition
+ if(!drainageFirst)
+ {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().Pcondition) continue;
+ cell->info().p()=bndCondValue[3];
+ cell->info().isWRes = false;
+ cell->info().isNWRes= true;
+ cell->info().saturation=0.0;
+ }
+ }
+ if(solver->debugOut) {cout<<"----initializeReservoirs----"<<endl;}
+}
+
+
#endif //TwoPhaseFLOW
=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp 2014-10-22 09:15:54 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp 2014-11-03 18:59:34 +0000
@@ -81,6 +81,10 @@
double bisection(CellHandle cell, int j, double a, double b);
double computeTriRadian(double a, double b, double c);
double computeDeltaForce(CellHandle cell,int j, double rC);
+ void initialization();
+ void computeSolidLine();
+ void initializeReservoirs();
+
//FIXME, needs to trigger initSolver() Somewhere, else changing flow.debug or other similar things after first calculation has no effect
//FIXME, I removed indexing cells from inside UnsatEngine (SoluteEngine shouldl be ok (?)) in order to get pressure computed, problem is they are not indexed at all if flow is not calculated
@@ -95,6 +99,8 @@
((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
((bool,initialWetting,true,,"Initial wetting saturated (=true) or non-wetting saturated (=false)"))
((bool, isPhaseTrapped,true,,"If True, both phases can be entrapped by the other, which would correspond to snap-off. If false, both phases are always connected to their reservoirs, thus no snap-off."))
+ ((bool, drainageFirst, true,,"If true, activate drainage first (initial saturated), then imbibition; if false, activate imbibition first (initial unsaturated), then drainage."))
+
,/*TwoPhaseFlowEngineT()*/,
,
@@ -104,6 +110,8 @@
.def("getCellSaturation",&TwoPhaseFlowEngine::cellSaturation,"get saturation of one pore")
.def("setCellSaturation",&TwoPhaseFlowEngine::setCellSaturation,"change saturation of one pore")
.def("computeOnePhaseFlow",&TwoPhaseFlowEngine::computeOnePhaseFlow,"compute pressure and fluxes in the W-phase")
+ .def("initialization",&TwoPhaseFlowEngine::initialization,"Initialize invasion setup. Build network, compute pore geometry info and initialize reservoir boundary conditions. ")
+
)
DECLARE_LOGGER;
};
=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp 2014-10-22 09:15:54 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp 2014-11-03 18:59:34 +0000
@@ -17,36 +17,39 @@
{
double totalCellVolume;
protected:
- void initialDrainage();
+// void initialization();
public :
- void initialReservoirs();///only used for determining first entry pressure
+// void initialReservoirs();///only used for determining first entry pressure
// void initializeCellIndex();
+// void initializeReservoirs();
void computeTotalPoresVolume();
double computeCellInterfacialArea(CellHandle cell, int j, double rC);
- void computeSolidLine();
+// void computeSolidLine();
void computeFacetPoreForcesWithCache(bool onlyCache=false);
void computeCapillaryForce() {computeFacetPoreForcesWithCache(false);}
- void drainage();
+ void invasion();
///functions can be shared by two modes
- void drainageSingleCell(CellHandle cell, double pressure);
+ void invasionSingleCell(CellHandle cell);
void updatePressure();
- double getMinEntryValue();
+ double getMinDrainagePc();
+ double getMaxImbibitionPc();
double getSaturation(bool isSideBoundaryIncluded=false);
double getSpecificInterfacialArea();
+// void updataPoreSaturation();
- void drainage1();
+ void invasion1();
void updateReservoirs1();
- void initialWResBound();
- void initialNWResBound();
+// void initialWResBound();
+// void initialNWResBound();
void WResRecursion(CellHandle cell);
void NWResRecursion(CellHandle cell);
- void checkWTrap(double pressure);
+ void checkTrap(double pressure);
- void drainage2();
+ void invasion2();
void updateReservoirs2();
void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
@@ -68,6 +71,8 @@
double getRMin(CellHandle cell, int j);
double getRMax(CellHandle cell, int j);
void checkRCompare();
+ void checkPoreRadius();
+ void printSomething();
virtual ~UnsaturatedEngine();
@@ -78,13 +83,15 @@
((bool, computeForceActivated, true,,"Activate capillary force computation. WARNING: turning off means capillary force is not computed at all, but the drainage can still work."))
((bool, isInvadeBoundary, true,,"Invasion side boundary condition. If True, pores of side boundary can be invaded; if False, the pore throats connecting side boundary are closed, those pores are excluded in saturation calculation."))
((int, windowsNo, 10,, "Number of genrated windows(or zoomed samples)."))
+// ((bool, drainageFirst, true,,"If true, activate drainage first (initial saturated), then imbibition; if false, activate imbibition first (initial unsaturated), then drainage."))
,,,
.def("saveVtk",&UnsaturatedEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
.def("savePhaseVtk",&UnsaturatedEngine::savePhaseVtk,(boost::python::arg("folder")="./phaseVtk"),"Save phases information in vtk format. W-phase=1, NW-phase=0, which can also be seen as saturation of local pore. Specify a folder name for output.")
- .def("getMinEntryValue",&UnsaturatedEngine::getMinEntryValue,"get the minimum NW entry pressure for the next invade step.")
+ .def("getMinDrainagePc",&UnsaturatedEngine::getMinDrainagePc,"Get the minimum entry capillary pressure for the next drainage step.")
+ .def("getMaxImbibitionPc",&UnsaturatedEngine::getMaxImbibitionPc,"Get the maximum entry capillary pressure for the next imbibition step.")
.def("getSaturation",&UnsaturatedEngine::getSaturation,(boost::python::arg("isSideBoundaryIncluded")),"Get saturation of entire packing. If isSideBoundaryIncluded=false (default), the pores of side boundary are excluded in saturation calculating; if isSideBoundaryIncluded=true (only in isInvadeBoundary=true drainage mode), the pores of side boundary are included in saturation calculating.")
.def("getSpecificInterfacialArea",&UnsaturatedEngine::getSpecificInterfacialArea,"get specific interfacial area (defined as the amount of fluid-fluid interfacial area per unit volume pf the porous medium).")
- .def("drainage",&UnsaturatedEngine::drainage,"Run the drainage invasion.")
+ .def("invasion",&UnsaturatedEngine::invasion,"Run the drainage invasion.")
.def("computeCapillaryForce",&UnsaturatedEngine::computeCapillaryForce,"Compute capillary force. ")
.def("checkCellsConnection",&UnsaturatedEngine::checkCellsConnection,"Check cell connections.")
@@ -95,11 +102,13 @@
.def("getInvadeDepth",&UnsaturatedEngine::getInvadeDepth,"Get NW-phase invasion depth. (the distance from NW-reservoir to front of NW-W interface.)")
.def("getSubdomainSaturation",&UnsaturatedEngine::getSubdomainSaturation,(boost::python::arg("pos"),boost::python::arg("radius")),"Get saturation of spherical subdomain defined by (pos, radius). The subdomain exclude boundary pores.")
- .def("initialDrainage",&UnsaturatedEngine::initialDrainage,"Initialize drainage setup. Build network, compute pore geometry info and initialize reservoir boundary conditions. ")
+// .def("initialization",&UnsaturatedEngine::initialization,"Initialize drainage setup. Build network, compute pore geometry info and initialize reservoir boundary conditions. ")
.def("checknoCache",&UnsaturatedEngine::checknoCache,"check noCache. (temp func.)")
.def("getWindowsSaturation",&UnsaturatedEngine::getWindowsSaturation,(boost::python::arg("windowsID"),boost::python::arg("isSideBoundaryIncluded")), "get saturation of subdomain with windowsID. If isSideBoundaryIncluded=false (default), the pores of side boundary are excluded in saturation calculating; if isSideBoundaryIncluded=true (only in isInvadeBoundary=true drainage mode), the pores of side boundary are included in saturation calculating.")
.def("initializeCellWindowsID",&UnsaturatedEngine::initializeCellWindowsID,"Initialize cell windows index. A temp function for comparison with experiments, will delete soon")
.def("checkRCompare",&UnsaturatedEngine::checkRCompare,"debug R RMin RMax.")
+ .def("checkPoreRadius",&UnsaturatedEngine::checkPoreRadius,"debug pore throat radius and pore body radius.")
+ .def("printSomething",&UnsaturatedEngine::printSomething,"print debug.")
)
DECLARE_LOGGER;
};
@@ -127,6 +136,7 @@
}
solver->noCache = true;
}*/
+/*
void UnsaturatedEngine::initialDrainage()
{
scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene
@@ -139,7 +149,22 @@
computeSolidLine();//save cell->info().solidLine[j][y]
initialReservoirs();
solver->noCache = true;
-}
+}*/
+/*
+void UnsaturatedEngine::initialization()
+{
+ scene = Omega::instance().getScene().get();//here define the pointer to Yade's scene
+ setPositionsBuffer(true);//copy sphere positions in a buffer...
+ buildTriangulation(0.0,*solver);//create a triangulation and initialize pressure in the elements (connecting with W-reservoir), everything will be contained in "solver"
+// initializeCellIndex();//initialize cell index
+ computePoreThroatRadius();//save pore throat radius before drainage
+ computePoreBodyRadius();//save pore body radius before imbibition
+ computeTotalPoresVolume();//save total volume of porous medium, considering different invading boundaries condition (isInvadeBoundary==True or False), aiming to calculate specific interfacial area.
+ computePoreBodyVolume();//save capillary volume of all cells, for fast calculating saturation
+ computeSolidLine();//save cell->info().solidLine[j][y]
+ initializeReservoirs();//initial pressure, reservoir flags and local pore saturation
+ solver->noCache = true;
+}*/
void UnsaturatedEngine::action()
{
@@ -190,7 +215,52 @@
if (cell->info().isFictious) continue;
totalCellVolume = totalCellVolume + std::abs( cell->info().volume() );}}
}
+/*
+void UnsaturatedEngine::initializeReservoirs()
+{
+ boundaryConditions(*solver);
+ solver->pressureChanged=true;
+ solver->reApplyBoundaryConditions();
+ ///keep boundingCells[2] as W-reservoir.
+ for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
+ (*it)->info().isWRes = true;
+ (*it)->info().isNWRes = false;
+ (*it)->info().saturation=1.0;
+ }
+ ///keep boundingCells[3] as NW-reservoir.
+ for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
+ (*it)->info().isNWRes = true;
+ (*it)->info().isWRes = false;
+ (*it)->info().saturation=0.0;
+ }
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ ///if we start from drainage
+ if(drainageFirst)
+ {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().Pcondition) continue;
+ cell->info().p()=bndCondValue[2];
+ cell->info().isWRes = true;
+ cell->info().isNWRes= false;
+ cell->info().saturation=1.0;
+ }
+ }
+ ///if we start from imbibition
+ if(!drainageFirst)
+ {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().Pcondition) continue;
+ cell->info().p()=bndCondValue[3];
+ cell->info().isWRes = false;
+ cell->info().isNWRes= true;
+ cell->info().saturation=0.0;
+ }
+ }
+ if(solver->debugOut) {cout<<"----initializeReservoirs----"<<endl;}
+}*/
+/*
void UnsaturatedEngine::initialReservoirs()
{
initialWResBound();
@@ -201,7 +271,7 @@
if (cell->info().p()==bndCondValue[2]) {cell->info().isWRes=true;cell->info().isNWRes=false;}
if (cell->info().p()==bndCondValue[3]) {cell->info().isNWRes=true;cell->info().isWRes=false;}
}
-}
+}*/
void UnsaturatedEngine::updatePressure()
{
@@ -214,12 +284,23 @@
if (cell->info().isWRes==true) {cell->info().p()=bndCondValue[2];}
if (cell->info().isNWRes==true) {cell->info().p()=bndCondValue[3];}
if (isPhaseTrapped) {
-// if ( (cell->info().isWRes==false)&&(cell->info().isNWRes==false) ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
if ( cell->info().isTrapW ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
+ if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
}
}
}
-
+/*
+void UnsaturatedEngine::updataPoreSaturation()
+{
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if( (cell->info().isWRes) || (cell->info().isTrapW) ) {cell->info().saturation=1.0;}
+ else if( (cell->info().isNWRes) || (cell->info().isTrapNW) ) {cell->info().saturation=0.0;}
+ else cerr<<"updataPoreSaturation Error! check pore flag."<<endl;
+ }
+}*/
+/*
///boundingCells[2] always connect W-reservoir.
void UnsaturatedEngine::initialWResBound()
{
@@ -233,14 +314,107 @@
for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
(*it)->info().isNWRes = true;
(*it)->info().isWRes = false;}
-}
-
-void UnsaturatedEngine::drainage()
-{
- if (isPhaseTrapped) drainage1();
- else drainage2();
-}
-
+}*/
+
+void UnsaturatedEngine::invasion()
+{
+ if (isPhaseTrapped) invasion1();
+ else invasion2();
+}
+
+///mode1 and mode2 can share the same invasionSingleCell(), invasionSingleCell() ONLY change neighbor pressure and neighbor saturation, independent of reservoirInfo.
+void UnsaturatedEngine::invasionSingleCell(CellHandle cell/*, double pressure*/)
+{
+ double localPressure=cell->info().p();
+ double localSaturation=cell->info().saturation;
+ for (int facet = 0; facet < 4; facet ++) {
+ CellHandle nCell = cell->neighbor(facet);
+ if (solver->T[solver->currentTes].Triangulation().is_infinite(nCell)) continue;
+ if (nCell->info().Pcondition) continue;//FIXME:defensive
+ if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+
+// if ( (nCell->info().saturation==localSaturation) && (nCell->info().p()==localPressure) ) continue;
+ if ( (nCell->info().saturation==localSaturation) && (nCell->info().p() != localPressure) && ((nCell->info().isTrapNW)||(nCell->info().isTrapW)) ) {
+ nCell->info().p() = localPressure;
+ if(solver->debugOut) {cerr<<"merge trapped phase"<<endl;}
+ invasionSingleCell(nCell);} ///here we merge trapped phase back to reservoir
+ else if ( (nCell->info().saturation>localSaturation) ) {
+ double nPcThroat=surfaceTension/cell->info().poreThroatRadius[facet];
+ double nPcBody=surfaceTension/nCell->info().poreBodyRadius;
+ if( (localPressure-nCell->info().p()>nPcThroat) && (localPressure-nCell->info().p()>nPcBody) ) {
+// if( localPressure-nCell->info().p()>nPcThroat ) {
+ nCell->info().p() = localPressure;
+ nCell->info().saturation=localSaturation;
+ if(solver->debugOut) {cerr<<"drainage"<<endl;}
+ invasionSingleCell(nCell);
+ }
+ }
+ else if ( (nCell->info().saturation<localSaturation) ) {
+// double nPcThroat=surfaceTension/nCell->info().poreThroatRadius[facet];
+ double nPcBody=surfaceTension/nCell->info().poreBodyRadius;
+ if( (nCell->info().p()-localPressure<nPcBody) && (nCell->info().p()-localPressure<surfaceTension/nCell->info().poreThroatRadius[0]) && (nCell->info().p()-localPressure<surfaceTension/nCell->info().poreThroatRadius[1]) && (nCell->info().p()-localPressure<surfaceTension/nCell->info().poreThroatRadius[2]) && (nCell->info().p()-localPressure<surfaceTension/nCell->info().poreThroatRadius[3]) ) {
+ nCell->info().p() = localPressure;
+ nCell->info().saturation=localSaturation;
+ if(solver->debugOut) {cerr<<"imbibition"<<endl;}
+ invasionSingleCell(nCell);
+ }
+ }
+ else continue;
+ }
+}
+///invasion mode 1: withTrap
+void UnsaturatedEngine::invasion1()
+{
+ if(solver->debugOut) {cout<<"----start drainage1----"<<endl;}
+
+ ///update Pw, Pn according to reservoirInfo.
+ updatePressure();
+ if(solver->debugOut) {cout<<"----invasion1.updatePressure----"<<endl;}
+
+// ///drainageSingleCell by Pressure difference, only change Pressure.
+// for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
+// invasionSingleCell((*it));
+// }
+// for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
+// invasionSingleCell((*it));
+// }
+// if(solver->debugOut) {cout<<"----invasion1.invasionSingleCell----"<<endl;}
+
+ ///drainageSingleCell by Pressure difference, only change Pressure.
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isWRes)
+ invasionSingleCell(cell);
+ }
+ ///drainageSingleCell by Pressure difference, only change Pressure.
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isNWRes)
+ invasionSingleCell(cell);
+ }
+
+ if(solver->debugOut) {cout<<"----invasion1.invasionSingleCell----"<<endl;}
+
+ ///update W, NW reservoirInfo according Pressure, trapped W-phase is marked by isWRes=False&&isNWRes=False.
+ updateReservoirs1();
+ if(solver->debugOut) {cout<<"----invasion1.update W, NW reservoirInfo----"<<endl;}
+
+ ///search new trapped W-phase/NW-phase, assign trapCapP for trapped W-phase
+ checkTrap(bndCondValue[3]-bndCondValue[2]);
+ if(solver->debugOut) {cout<<"----invasion1.checkWTrap----"<<endl;}
+
+ ///update trapped W-phase/NW-phase Pressure
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if ( cell->info().isTrapW ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
+ if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
+ }
+ if(solver->debugOut) {cout<<"----invasion1.update trapped W-phase/NW-phase Pressure----"<<endl;}
+
+// ///update local pore saturation
+// updataPoreSaturation();
+// if(solver->debugOut) {cout<<"----drainage1.update local pore saturation----"<<endl;}
+}
+/*
///mode1 and mode2 can share the same drainageSingleCell(), drainageSingleCell() ONLY change pressure, independent of reservoirInfo.
void UnsaturatedEngine::drainageSingleCell(CellHandle cell, double pressure)
{
@@ -254,8 +428,8 @@
if (pressure-nCell->info().p() > nCellP) {
nCell->info().p() = pressure;
drainageSingleCell(nCell, pressure);}}}
-}
-
+}*/
+/*
///drainage mode 1: withTrap
void UnsaturatedEngine::drainage1()
{
@@ -283,26 +457,28 @@
if(solver->debugOut) {cout<<"----drainage1.checkWTrap----"<<endl;}
///update trapped W-phase Pressure
- FiniteCellsIterator ncellEnd = tri.finite_cells_end();
- for ( FiniteCellsIterator ncell = tri.finite_cells_begin(); ncell != ncellEnd; ncell++ ) {
-// if( (ncell->info().isWRes) || (ncell->info().isNWRes) ) continue;
- if(ncell->info().isTrapW)
- ncell->info().p() = bndCondValue[3] - ncell->info().trapCapP;
+ FiniteCellsIterator nCellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator nCell = tri.finite_cells_begin(); nCell != nCellEnd; nCell++ ) {
+ if(nCell->info().isTrapW)
+ nCell->info().p() = bndCondValue[3] - nCell->info().trapCapP;
}
if(solver->debugOut) {cout<<"----drainage1.update trapped W-phase Pressure----"<<endl;}
-}
+ ///update local pore saturation
+ updataPoreSaturation();
+ if(solver->debugOut) {cout<<"----drainage1.update local pore saturation----"<<endl;}
+}*/
-///search trapped W-phase, define trapCapP=Pn-Pw. assign isTrapW info.
-void UnsaturatedEngine::checkWTrap(double pressure)
+///search trapped W-phase or NW-phase, define trapCapP=Pn-Pw. assign isTrapW/isTrapNW info.
+void UnsaturatedEngine::checkTrap(double pressure)
{
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if( (cell->info().isWRes) || (cell->info().isNWRes) || (cell->info().isTrapW) ) continue;
-// if(cell->info().p()!= bndCondValue[2]) continue;
+ if( (cell->info().isWRes) || (cell->info().isNWRes) || (cell->info().isTrapW) || (cell->info().isTrapNW) ) continue;
cell->info().trapCapP=pressure;
- cell->info().isTrapW=true;
+ if(cell->info().saturation==1.0) cell->info().isTrapW=true;
+ if(cell->info().saturation==0.0) cell->info().isTrapNW=true;
}
}
@@ -311,12 +487,13 @@
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().Pcondition) continue;
cell->info().isWRes = false;
cell->info().isNWRes = false;
}
-
- initialWResBound();
- initialNWResBound();
+//
+// initialWResBound();
+// initialNWResBound();
for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
if ((*it)==NULL) continue;
@@ -340,6 +517,8 @@
if (nCell->info().isWRes==true) continue;
nCell->info().isWRes = true;
nCell->info().isNWRes = false;
+ nCell->info().isTrapW = false;
+ nCell->info().trapCapP=0.0;
WResRecursion(nCell);
}
}
@@ -355,24 +534,40 @@
if (nCell->info().isNWRes==true) continue;
nCell->info().isNWRes = true;
nCell->info().isWRes = false;
+ nCell->info().isTrapNW = false;
+ nCell->info().trapCapP=0.0;
NWResRecursion(nCell);
}
}
-///drainage mode 2: withoutTrap
-void UnsaturatedEngine::drainage2()
+///invasion mode 2: withoutTrap
+void UnsaturatedEngine::invasion2()
{
+ if(solver->debugOut) {cout<<"----start invasion2----"<<endl;}
+
///update Pw, Pn according to reservoirInfo.
updatePressure();
+ if(solver->debugOut) {cout<<"----invasion2.updatePressure----"<<endl;}
+
///drainageSingleCell by Pressure difference, only change Pressure.
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if(cell->info().p() == bndCondValue[3])
- drainageSingleCell(cell,cell->info().p());
- }
+ if(cell->info().isWRes)
+ invasionSingleCell(cell);
+ }
+ ///drainageSingleCell by Pressure difference, only change Pressure.
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isNWRes)
+ invasionSingleCell(cell);
+ }
+
+ if(solver->debugOut) {cout<<"----invasion2.invasionSingleCell----"<<endl;}
+
///update W, NW reservoirInfo according Pressure
updateReservoirs2();
+ if(solver->debugOut) {cout<<"----drainage2.update W, NW reservoirInfo----"<<endl;}
+
}
void UnsaturatedEngine::updateReservoirs2()
@@ -386,7 +581,7 @@
}
}
-double UnsaturatedEngine::getMinEntryValue()
+double UnsaturatedEngine::getMinDrainagePc()
{
double nextEntry = 1e50;
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
@@ -394,11 +589,13 @@
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
if (cell->info().isNWRes == true) {
for (int facet=0; facet<4; facet ++) {
- if (tri.is_infinite(cell->neighbor(facet))) continue;
- if (cell->neighbor(facet)->info().Pcondition) continue;
- if ( (cell->neighbor(facet)->info().isFictious) && (!isInvadeBoundary) ) continue;
- if ( cell->neighbor(facet)->info().isWRes == true ) {
- double nCellP = surfaceTension/cell->info().poreThroatRadius[facet];
+ CellHandle nCell = cell->neighbor(facet);
+ if (tri.is_infinite(nCell)) continue;
+ if (nCell->info().Pcondition) continue;
+ if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if ( nCell->info().isWRes == true ) {
+ double nCellP = std::max( (surfaceTension/cell->info().poreThroatRadius[facet]),(surfaceTension/nCell->info().poreBodyRadius) );
+// double nCellP = surfaceTension/cell->info().poreThroatRadius[facet];
nextEntry = std::min(nextEntry,nCellP);}}}}
if (nextEntry==1e50) {
@@ -408,6 +605,31 @@
else return nextEntry;
}
+double UnsaturatedEngine::getMaxImbibitionPc()
+{
+ double nextEntry = -1e50;
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().isWRes == true) {
+ for (int facet=0; facet<4; facet ++) {
+ CellHandle nCell = cell->neighbor(facet);
+ if (tri.is_infinite(nCell)) continue;
+ if (nCell->info().Pcondition) continue;
+ if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if ( nCell->info().isNWRes == true ) {
+// double nCellP = std::min( (surfaceTension/nCell->info().poreBodyRadius), (surfaceTension/cell->info().poreThroatRadius[facet]));
+ double nCellP = std::min( (surfaceTension/nCell->info().poreBodyRadius), std::min( std::min( std::min( (surfaceTension/nCell->info().poreThroatRadius[0]),(surfaceTension/nCell->info().poreThroatRadius[1]) ), (surfaceTension/nCell->info().poreThroatRadius[2]) ),(surfaceTension/nCell->info().poreThroatRadius[3])) );
+ nextEntry = std::max(nextEntry,nCellP);}}}}
+
+ if (nextEntry==-1e50) {
+ cout << "End imbibition !" << endl;
+ return nextEntry=0;
+ }
+ else return nextEntry;
+}
+
+/*
double UnsaturatedEngine::getSaturation(bool isSideBoundaryIncluded)
{
if( (!isInvadeBoundary) && (isSideBoundaryIncluded)) cerr<<"In isInvadeBoundary=false drainage, isSideBoundaryIncluded can't set true."<<endl;
@@ -426,6 +648,25 @@
}
double saturation = 1 - nwVolume/poresVolume;
return saturation;
+}*/
+
+double UnsaturatedEngine::getSaturation(bool isSideBoundaryIncluded)
+{
+ if( (!isInvadeBoundary) && (isSideBoundaryIncluded)) cerr<<"In isInvadeBoundary=false drainage, isSideBoundaryIncluded can't set true."<<endl;
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ double poresVolume = 0.0; //total pores volume
+ double wVolume = 0.0; //NW-phase volume
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().Pcondition) continue;
+ if ( (cell->info().isFictious) && (!isSideBoundaryIncluded) ) continue;
+ poresVolume = poresVolume + cell->info().poreBodyVolume;
+ if (cell->info().saturation>0.0) {
+ wVolume = wVolume + cell->info().poreBodyVolume * cell->info().saturation;
+ }
+ }
+ return wVolume/poresVolume;
}
double UnsaturatedEngine::getSpecificInterfacialArea()
@@ -546,15 +787,49 @@
}
vtkfile.end_cells();
- vtkfile.begin_data("Phase",CELL_DATA,SCALARS,FLOAT);
+ vtkfile.begin_data("Saturation",CELL_DATA,SCALARS,FLOAT);
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
bool isDrawable = cell->info().isReal() && cell->vertex(0)->info().isReal() && cell->vertex(1)->info().isReal() && cell->vertex(2)->info().isReal() && cell->vertex(3)->info().isReal();
- if (isDrawable){vtkfile.write_data(!cell->info().isNWRes);}
+ if (isDrawable){vtkfile.write_data(cell->info().saturation);}
}
vtkfile.end_data();
}
// ------------------for checking----------
+void UnsaturatedEngine::checkPoreRadius()
+{
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ ofstream file;
+ file.open("checkPoreRadius.txt");
+ file<<"cellID "<<"poreThroatRadius[0] "<<"poreThroatRadius[1] "<<"poreThroatRadius[2] "<<"poreThroatRadius[3] "<<"poreBodyRadius"<<endl;
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ file << cell->info().index << " " <<cell->info().poreThroatRadius[0]<< " "<<cell->info().poreThroatRadius[1]<< " "<<cell->info().poreThroatRadius[2]<< " "<<cell->info().poreThroatRadius[3]<< " "<<cell->info().poreBodyRadius<< endl;
+ }
+ file.close();
+ file.open("checkPoreThroatRadiusCompareBodyRadius.txt");
+ file<<"cellID "<<"poreThroatRadius[0] "<<"poreThroatRadius[1] "<<"poreThroatRadius[2] "<<"poreThroatRadius[3] "<<"poreBodyRadius"<<endl;
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if( (cell->info().poreBodyRadius<cell->info().poreThroatRadius[0]) || (cell->info().poreBodyRadius<cell->info().poreThroatRadius[1]) || (cell->info().poreBodyRadius<cell->info().poreThroatRadius[2]) || (cell->info().poreBodyRadius<cell->info().poreThroatRadius[3]) )
+ file << cell->info().index << " " <<cell->info().poreThroatRadius[0]<< " "<<cell->info().poreThroatRadius[1]<< " "<<cell->info().poreThroatRadius[2]<< " "<<cell->info().poreThroatRadius[3]<< " "<<cell->info().poreBodyRadius<< endl;
+ }
+ file.close();
+ file.open("checkPoreThroatCircleRadiusCompareBodyRadius.txt");
+ file<<"cellID "<<"poreThroatRadius[0] "<<"poreThroatRadius[1] "<<"poreThroatRadius[2] "<<"poreThroatRadius[3] "<<"poreBodyRadius"<<endl;
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if( (cell->info().poreBodyRadius<solver->computeEffectiveRadius(cell,0)) || (cell->info().poreBodyRadius<solver->computeEffectiveRadius(cell,1)) || (cell->info().poreBodyRadius<solver->computeEffectiveRadius(cell,2)) || (cell->info().poreBodyRadius<solver->computeEffectiveRadius(cell,3)) )
+ file << cell->info().index << " " <<solver->computeEffectiveRadius(cell,0)<< " "<<solver->computeEffectiveRadius(cell,1)<< " "<<solver->computeEffectiveRadius(cell,2)<< " "<<solver->computeEffectiveRadius(cell,3)<< " "<<cell->info().poreBodyRadius<< endl;
+ }
+ file.close();
+ file.open("checkPoreBodyRadiusNegative.txt");
+ file<<"cellID "<<"poreThroatRadius[0] "<<"poreThroatRadius[1] "<<"poreThroatRadius[2] "<<"poreThroatRadius[3] "<<"poreBodyRadius"<<endl;
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if( (cell->info().poreBodyRadius<0) )
+ file << cell->info().index << " " <<solver->computeEffectiveRadius(cell,0)<< " "<<solver->computeEffectiveRadius(cell,1)<< " "<<solver->computeEffectiveRadius(cell,2)<< " "<<solver->computeEffectiveRadius(cell,3)<< " "<<cell->info().poreBodyRadius<< endl;
+ }
+ file.close();
+
+}
double UnsaturatedEngine::getRMin(CellHandle cell, int j)
{
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
@@ -665,15 +940,23 @@
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
ofstream file;
file.open("entryCapillaryPressure.txt");
- file << "#List of Connections \n";
- file << "cell" << " " << "neighborCell" << " " << "entryCapillaryPressure" << endl;
FiniteCellsIterator cellEnd = tri.finite_cells_end();
for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
for (int facet=0; facet<4; facet++) {
- file << cell->info().index << " " <<cell->neighbor(facet)->info().index << " " << surfaceTension/cell->info().poreThroatRadius[facet] << endl;
- }
- }
- file.close();
+ file <<surfaceTension/cell->info().poreThroatRadius[facet] << endl;
+ }
+ }
+ file.close();
+
+ file.open("entryCapillaryPressureNoBound.txt");
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if (cell->info().isFictious) continue;
+ for (int facet=0; facet<4; facet++) {
+ file <<surfaceTension/cell->info().poreThroatRadius[facet] << endl;
+ }
+ }
+ file.close();
+
}
void UnsaturatedEngine::checkLatticeNodeY(double y)
@@ -755,6 +1038,17 @@
}
}
+void UnsaturatedEngine::printSomething()
+{
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if ((cell->info().isFictious==true)&&(cell->info().Pcondition==false)) {
+cerr<<cell->info().saturation<<endl; }
+ }
+
+}
+
//----------temp functions for comparison with experiment-----------------------
void UnsaturatedEngine::initializeCellWindowsID()
{
@@ -772,7 +1066,7 @@
{
if( (!isInvadeBoundary) && (isSideBoundaryIncluded)) cerr<<"In isInvadeBoundary=false drainage, isSideBoundaryIncluded can't set true."<<endl;
double poresVolume = 0.0; //total pores volume
- double nwVolume = 0.0; //NW-phase volume
+ double wVolume = 0.0; //W-phase volume
RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = tri.finite_cells_end();
@@ -781,12 +1075,11 @@
if ( (cell->info().isFictious) && (!isSideBoundaryIncluded) ) continue;
if (cell->info().windowsID != i) continue;
poresVolume = poresVolume + cell->info().poreBodyVolume;
- if (cell->info().p()==bndCondValue[3]) {
- nwVolume = nwVolume + cell->info().poreBodyVolume;
+ if (cell->info().saturation>0.0) {
+ wVolume = wVolume + cell->info().poreBodyVolume * cell->info().saturation;
}
}
- double saturation = 1 - nwVolume/poresVolume;
- return saturation;
+ return wVolume/poresVolume;
}
double UnsaturatedEngine::getInvadeDepth()
@@ -807,9 +1100,8 @@
double UnsaturatedEngine::getSubdomainSaturation(Vector3r pos, double radius)
{
- double sw=0.0;
double poresVolume=0.0;
- double nwVolume=0.0;
+ double wVolume=0.0;
RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
FiniteCellsIterator cellEnd = Tri.finite_cells_end();
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
@@ -817,24 +1109,21 @@
double dist=(pos-cellPos).norm();
if(dist>radius) continue;
if(cell->info().isFictious) {
- sw=-1;
- break;
+ cerr<<"The radius of subdomain is too large, or the center of subdomain is out of packing. Please reset subdomain again."<<endl;
+ return -1;
}
poresVolume=poresVolume+cell->info().poreBodyVolume;
- if(cell->info().isNWRes) {
- nwVolume=nwVolume+cell->info().poreBodyVolume;
+ if(cell->info().saturation>0.0) {
+ wVolume=wVolume+cell->info().poreBodyVolume * cell->info().saturation;
}
}
- if(sw==-1) {
- cerr<<"The radius of subdomain is too large, or the center of subdomain is out of packing. Please reset subdomain again."<<endl;
- }
- else sw=1-nwVolume/poresVolume;
- return sw;
+ return wVolume/poresVolume;
}
//--------------end of comparison with experiment----------------------------
///compute forces
+/*
void UnsaturatedEngine::computeSolidLine()
{
RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
@@ -845,7 +1134,7 @@
}
}
if(solver->debugOut) {cout<<"----computeSolidLine-----."<<endl;}
-}
+}*/
void UnsaturatedEngine::computeFacetPoreForcesWithCache(bool onlyCache)
{
