yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12932
[Branch ~yade-pkg/yade/git-trunk] Rev 3975: clean UnsaturatedEngine, move all functions to TwoPhaseFlowEngine.
------------------------------------------------------------
revno: 3975
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2016-11-22 18:29:15 +0100
message:
clean UnsaturatedEngine, move all functions to TwoPhaseFlowEngine.
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 2016-11-02 10:40:55 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp 2016-11-22 17:29:15 +0000
@@ -1,4 +1,3 @@
-
/*************************************************************************
* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2013 by T. Sweijen (T.sweijen@xxxxx) *
@@ -43,12 +42,6 @@
solver->noCache = true;
}
-
-
-
-
-
-
void TwoPhaseFlowEngine::computePoreBodyVolume()
{
initializeVolumes(*solver);
@@ -59,8 +52,6 @@
}
}
-
-
void TwoPhaseFlowEngine::computePoreThroatRadiusMethod2()
{
//Calculate the porethroat radii of the inscribed sphere in each pore-body.
@@ -85,8 +76,6 @@
}
}
-
-
void TwoPhaseFlowEngine::computePoreBodyRadius()
{
@@ -435,8 +424,6 @@
}
}
-
-
void TwoPhaseFlowEngine::computeSolidLine()
{
RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
@@ -713,5 +700,420 @@
}
}
+void TwoPhaseFlowEngine::updatePressure()
+{
+ boundaryConditions(*solver);
+ solver->pressureChanged=true;
+ solver->reApplyBoundaryConditions();
+ 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) {cell->info().p()=bndCondValue[2];}
+ if (cell->info().isNWRes==true) {cell->info().p()=bndCondValue[3];}
+ if (isPhaseTrapped) {
+ if ( cell->info().isTrapW ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
+ if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
+ //check cell reservoir info.
+ if ( !cell->info().isWRes && !cell->info().isNWRes && !cell->info().isTrapW && !cell->info().isTrapNW ) {cerr<<"ERROR! NOT FIND Cell Info!";}
+// {cell->info().p()=bndCondValue[2]; if (isInvadeBoundary) cerr<<"Something wrong in updatePressure.(isInvadeBoundary)";}
+ }
+ }
+}
+
+void TwoPhaseFlowEngine::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 TwoPhaseFlowEngine::invasionSingleCell(CellHandle cell)
+{
+ 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 (cell->info().poreThroatRadius[facet]<0) 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) ) {
+ nCell->info().p() = localPressure;
+ nCell->info().saturation=localSaturation;
+ nCell->info().hasInterface=false;
+ if(solver->debugOut) {cerr<<"drainage"<<endl;}
+ if (recursiveInvasion) invasionSingleCell(nCell);
+ }
+////FIXME:Introduce cell.hasInterface
+// else if( (localPressure-nCell->info().p()>nPcThroat) && (localPressure-nCell->info().p()<nPcBody) && (cell->info().hasInterface==false) && (nCell->info().hasInterface==false) ) {
+// if(solver->debugOut) {cerr<<"invasion paused into pore interface "<<endl;}
+// nCell->info().hasInterface=true;
+// }
+// else continue;
+ }
+ else if ( (nCell->info().saturation<localSaturation) ) {
+ double nPcThroat=surfaceTension/cell->info().poreThroatRadius[facet];
+ double nPcBody=surfaceTension/nCell->info().poreBodyRadius;
+ if( (nCell->info().p()-localPressure<nPcBody) && (nCell->info().p()-localPressure<nPcThroat) ) {
+ nCell->info().p() = localPressure;
+ nCell->info().saturation=localSaturation;
+ if(solver->debugOut) {cerr<<"imbibition"<<endl;}
+ if (recursiveInvasion) invasionSingleCell(nCell);
+ }
+//// FIXME:Introduce cell.hasInterface
+// else if ( (nCell->info().p()-localPressure<nPcBody) && (nCell->info().p()-localPressure>nPcThroat) /*&& (cell->info().hasInterface==false) && (nCell->info().hasInterface==false)*/ ) {
+// nCell->info().p() = localPressure;
+// nCell->info().saturation=localSaturation;
+// if(solver->debugOut) {cerr<<"imbibition paused pore interface"<<endl;}
+// nCell->info().hasInterface=true;
+// }
+// else continue;
+ }
+ else continue;
+ }
+}
+///invasion mode 1: withTrap
+void TwoPhaseFlowEngine::invasion1()
+{
+ if(solver->debugOut) {cout<<"----start invasion1----"<<endl;}
+
+ ///update Pw, Pn according to reservoirInfo.
+ updatePressure();
+ if(solver->debugOut) {cout<<"----invasion1.updatePressure----"<<endl;}
+
+ ///invasionSingleCell by Pressure difference, change Pressure and Saturation.
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ if(isDrainageActivated) {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isNWRes)
+ invasionSingleCell(cell);
+ }
+ }
+ if(isImbibitionActivated) {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isWRes)
+ invasionSingleCell(cell);
+ }
+ }
+ if(solver->debugOut) {cout<<"----invasion1.invasionSingleCell----"<<endl;}
+
+ ///update W, NW reservoirInfo according to cell->info().saturation
+ updateReservoirs1();
+ if(solver->debugOut) {cout<<"----invasion1.update W, NW reservoirInfo----"<<endl;}
+
+ ///search new trapped W-phase/NW-phase, assign trapCapP, isTrapW/isTrapNW flag for new trapped phases. But at this moment, the new trapped W/NW cells.P= W/NW-Res.P. They will be updated in next updatePressure() func.
+ checkTrap(bndCondValue[3]-bndCondValue[2]);
+ if(solver->debugOut) {cout<<"----invasion1.checkWTrap----"<<endl;}
+
+ ///update trapped W-phase/NW-phase Pressure //FIXME: is this necessary?
+ 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;}
+
+ if(isCellLabelActivated) updateCellLabel();
+ if(solver->debugOut) {cout<<"----update cell labels----"<<endl;}
+}
+
+///search trapped W-phase or NW-phase, define trapCapP=Pn-Pw. assign isTrapW/isTrapNW info.
+void TwoPhaseFlowEngine::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().isFictious) && (!cell->info().Pcondition) && (!isInvadeBoundary) ) continue;
+ if( (cell->info().isWRes) || (cell->info().isNWRes) || (cell->info().isTrapW) || (cell->info().isTrapNW) ) continue;
+ cell->info().trapCapP=pressure;
+ if(cell->info().saturation==1.0) cell->info().isTrapW=true;
+ if(cell->info().saturation==0.0) cell->info().isTrapNW=true;
+ }
+}
+
+void TwoPhaseFlowEngine::updateReservoirs1()
+{
+ 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;
+ }
+
+ for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
+ if ((*it)==NULL) continue;
+ WResRecursion(*it);
+ }
+
+ for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
+ if ((*it)==NULL) continue;
+ NWResRecursion(*it);
+ }
+}
+
+void TwoPhaseFlowEngine::WResRecursion(CellHandle cell)
+{
+ 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;
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if (nCell->info().saturation != 1.0) continue;
+ 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);
+ }
+}
+
+void TwoPhaseFlowEngine::NWResRecursion(CellHandle cell)
+{
+ 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;
+// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
+ if (nCell->info().saturation != 0.0) continue;
+ 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);
+ }
+}
+
+///invasion mode 2: withoutTrap
+void TwoPhaseFlowEngine::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, change Pressure and Saturation.
+ RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
+ FiniteCellsIterator cellEnd = tri.finite_cells_end();
+ if(isDrainageActivated) {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isNWRes)
+ invasionSingleCell(cell);
+ }
+ }
+ ///drainageSingleCell by Pressure difference, change Pressure and Saturation.
+ if(isImbibitionActivated) {
+ for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
+ if(cell->info().isWRes)
+ invasionSingleCell(cell);
+ }
+ }
+
+ if(solver->debugOut) {cout<<"----invasion2.invasionSingleCell----"<<endl;}
+
+ ///update W, NW reservoirInfo according to Pressure
+ updateReservoirs2();
+ if(solver->debugOut) {cout<<"----drainage2.update W, NW reservoirInfo----"<<endl;}
+
+}
+
+void TwoPhaseFlowEngine::updateReservoirs2()
+{
+ 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[2]) {cell->info().isWRes=true; cell->info().isNWRes=false;}
+ else if (cell->info().p()==bndCondValue[3]) {cell->info().isNWRes=true; cell->info().isWRes=false;}
+ else {cerr<<"drainage mode2: updateReservoir Error!"<<endl;}
+ }
+}
+
+double TwoPhaseFlowEngine::getMinDrainagePc()
+{
+ 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().isNWRes == 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().isWRes == true && cell->info().poreThroatRadius[facet]>0) {
+ 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) {
+ cout << "End drainage !" << endl;
+ return nextEntry=0;
+ }
+ else return nextEntry;
+}
+
+double TwoPhaseFlowEngine::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 && cell->info().poreThroatRadius[facet]>0) {
+ double nCellP = std::min( (surfaceTension/nCell->info().poreBodyRadius), (surfaceTension/cell->info().poreThroatRadius[facet]));
+ nextEntry = std::max(nextEntry,nCellP);}}}}
+
+ if (nextEntry==-1e50) {
+ cout << "End imbibition !" << endl;
+ return nextEntry=0;
+ }
+ else return nextEntry;
+}
+
+double TwoPhaseFlowEngine::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;
+}
+
+///compute forces
+void TwoPhaseFlowEngine::computeFacetPoreForcesWithCache(bool onlyCache)
+{
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ CVector nullVect(0,0,0);
+ //reset forces
+ if (!onlyCache) for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
+// #ifdef parallel_forces
+// if (solver->noCache) {
+// solver->perVertexUnitForce.clear(); solver->perVertexPressure.clear();
+// solver->perVertexUnitForce.resize(solver->T[solver->currentTes].maxId+1);
+// solver->perVertexPressure.resize(solver->T[solver->currentTes].maxId+1);}
+// #endif
+// CellHandle neighbourCell;
+// VertexHandle mirrorVertex;
+ CVector tempVect;
+ //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
+ if (solver->noCache) {//WARNING:all currentTes must be solver->T[solver->currentTes], should NOT be solver->T[currentTes]
+ for (FlowSolver::VCellIterator cellIt=solver->T[solver->currentTes].cellHandles.begin(); cellIt!=solver->T[solver->currentTes].cellHandles.end(); cellIt++) {
+ CellHandle& cell = *cellIt;
+ //reset cache
+ for (int k=0; k<4; k++) cell->info().unitForceVectors[k]=nullVect;
+
+ for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
+ const CVector& Surfk = cell->info().facetSurfaces[j];
+ //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
+ //The ratio void surface / facet surface
+ Real area = sqrt(Surfk.squared_length());
+ if (area<=0) cerr <<"AREA <= 0!! AREA="<<area<<endl;
+ CVector facetNormal = Surfk/area;
+ const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
+ CVector fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
+ /// handle fictious vertex since we can get the projected surface easily here
+ if (cell->vertex(j)->info().isFictious) {
+ Real projSurf=std::abs(Surfk[solver->boundary(cell->vertex(j)->info().id()).coordinate]);
+ tempVect=-projSurf*solver->boundary(cell->vertex(j)->info().id()).normal;
+ cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p();
+ //define the cached value for later use with cache*p
+ cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
+ }
+ /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
+ CVector facetUnitForce = -fluidSurfk*cell->info().solidLine[j][3];
+ CVector facetForce = cell->info().p()*facetUnitForce;
+
+ for (int y=0; y<3; y++) {
+ cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + facetForce*cell->info().solidLine[j][y];
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+facetUnitForce*cell->info().solidLine[j][y];
+ //uncomment to get total force / comment to get only pore tension forces
+ if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
+ cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces -facetNormal*cell->info().p()*crossSections[j][y];
+ //add to cached value
+ cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
+ }
+ }
+// #ifdef parallel_forces
+// solver->perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j]));
+// solver->perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p()));
+// #endif
+ }
+ }
+ solver->noCache=false;//cache should always be defined after execution of this function
+ }
+ if (onlyCache) return;
+
+// else {//use cached values when triangulation doesn't change
+// #ifndef parallel_forces
+ for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); cell++) {
+ for (int yy=0; yy<4; yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();
+ }
+
+// #else
+// #pragma omp parallel for num_threads(ompThreads)
+// for (int vn=0; vn<= solver->T[solver->currentTes].maxId; vn++) {
+// VertexHandle& v = solver->T[solver->currentTes].vertexHandles[vn];
+// const int& id = v->info().id();
+// CVector tf (0,0,0);
+// int k=0;
+// for (vector<const Real*>::iterator c = solver->perVertexPressure[id].begin(); c != solver->perVertexPressure[id].end(); c++)
+// tf = tf + (*(solver->perVertexUnitForce[id][k++]))*(**c);
+// v->info().forces = tf;
+// }
+// #endif
+// }
+ if (solver->debugOut) {
+ CVector totalForce = nullVect;
+ for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
+ if (!v->info().isFictious) totalForce = totalForce + v->info().forces;
+ else if (solver->boundary(v->info().id()).flowCondition==1) totalForce = totalForce + v->info().forces;
+ }
+ cout << "totalForce = "<< totalForce << endl;
+ }
+}
+
+bool TwoPhaseFlowEngine::detectBridge(RTriangulation::Finite_edges_iterator& edge)
+{
+ bool dryBridgeExist=true;
+ const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ RTriangulation::Cell_circulator cell1 = Tri.incident_cells(*edge);
+ RTriangulation::Cell_circulator cell0 = cell1++;
+ if(cell0->info().saturation==1) {dryBridgeExist=false; return dryBridgeExist;}
+ else {
+ while (cell1!=cell0) {
+ if(cell1->info().saturation==1) {dryBridgeExist=false;break;}
+ else cell1++;}
+ return dryBridgeExist;
+ }
+}
+
+
#endif //TwoPhaseFLOW
=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp 2016-11-02 10:47:10 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp 2016-11-22 17:29:15 +0000
@@ -1,4 +1,3 @@
-
/*************************************************************************
* Copyright (C) 2014 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2013 by Thomas. Sweijen <T.sweijen@xxxxx> *
@@ -67,38 +66,75 @@
//if we overload action() like this, this engine is doing nothing in a standard timestep, it can still have useful functions
virtual void action() {};
- //If a new function is specific to the derived engine, put it here, else go to the base TemplateFlowEngine
- //if it is useful for everyone
- void computePoreBodyVolume();
- void computePoreBodyRadius();
-// void computePoreThroatCircleRadius();
- double computePoreSatAtInterface(int ID);
- void computePoreCapillaryPressure(CellHandle cell);
- void savePhaseVtk(const char* folder);
-// void computePoreThroatRadius();
-// void computePoreThroatRadiusTricky();//set the radius of pore throat between side pores negative.
-
+ //If a new function is specific to the derived engine, put it here, else go to the base TemplateFlowEngine if it is useful for everyone
+
+ //compute pore body radius
+ void computePoreBodyVolume();
+ void computePoreBodyRadius();
+ void computeSolidLine();
+
+ //compute entry pore throat radius (drainage)
+ void computePoreThroatRadiusMethod1();//MS-P method
+ void computePoreThroatRadiusTrickyMethod1();//set the radius of pore throat between side pores negative.
double computeEffPoreThroatRadius(CellHandle cell, int j);
double computeEffPoreThroatRadiusFine(CellHandle cell, int j);
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 computePoreThroatRadiusMethod2();//radius of the inscribed circle
+ void computePoreThroatRadiusMethod3();//radius of area-equivalent circle
+
+ ///begin of invasion (mainly drainage) model
void initialization();
- void computeSolidLine();
void initializeReservoirs();
-
- void computePoreThroatRadiusMethod1();
- void computePoreThroatRadiusTrickyMethod1();//set the radius of pore throat between side pores negative.
- void computePoreThroatRadiusMethod2();
- void computePoreThroatRadiusMethod3();
- void savePoreNetwork();
-
+
+ void invasion();//functions can be shared by two modes
+ void invasionSingleCell(CellHandle cell);
+ void updatePressure();
+ double getMinDrainagePc();
+ double getMaxImbibitionPc();
+ double getSaturation(bool isSideBoundaryIncluded=false);
+
+ void invasion1();//with-trap
+ void updateReservoirs1();
+ void WResRecursion(CellHandle cell);
+ void NWResRecursion(CellHandle cell);
+ void checkTrap(double pressure);
void updateReservoirLabel();
void updateCellLabel();
void updateSingleCellLabelRecursion(CellHandle cell, int label);
int getMaxCellLabel();
-
+ void invasion2();//without-trap
+ void updateReservoirs2();
+ ///end of invasion model
+
+ //compute forces
+ void computeFacetPoreForcesWithCache(bool onlyCache=false);
+ void computeCapillaryForce() {computeFacetPoreForcesWithCache(false);}
+
+ //combine with pendular model
+ boost::python::list getPotentialPendularSpheresPair() {
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ boost::python::list bridgeIds;
+ FiniteEdgesIterator ed_it = Tri.finite_edges_begin();
+ for ( ; ed_it!=Tri.finite_edges_end(); ed_it++ ) {
+ if (detectBridge(ed_it)==true) {
+ const VertexInfo& vi1=(ed_it->first)->vertex(ed_it->second)->info();
+ const VertexInfo& vi2=(ed_it->first)->vertex(ed_it->third)->info();
+ const int& id1 = vi1.id();
+ const int& id2 = vi2.id();
+ bridgeIds.append(boost::python::make_tuple(id1,id2));}}
+ return bridgeIds;
+ }
+ bool detectBridge(RTriangulation::Finite_edges_iterator& edge);
+
+ //post-processing
+ void savePoreNetwork();
+ void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
+ void savePhaseVtk(const char* folder);
+
boost::python::list cellporeThroatRadius(unsigned int id){ // Temporary function to allow for simulations in Python, can be easily accessed in c++
boost::python::list ids;
if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
@@ -112,7 +148,10 @@
for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->neighbor(i)->info().id);
return ids;
}
-
+
+ //TODO
+ double computePoreSatAtInterface(int ID);
+ void computePoreCapillaryPressure(CellHandle cell);
//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
@@ -143,7 +182,11 @@
((int,entryPressureMethod,1,,"integer to define the method used to determine the pore throat radii and the according entry pressures. 1)radius of entry pore throat based on MS-P method; 2) radius of the inscribed circle; 3) radius of the circle with equivalent surface area of the pore throat."))
((double,partiallySaturatedPores,false,,"Include partially saturated pores or not?"))
((bool, isCellLabelActivated, false,, "Activate cell labels for marking disconnected wetting clusters. NW-reservoir label 0; W-reservoir label 1; disconnected W-clusters label from 2. "))
+ ((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, isDrainageActivated, true,, "Activates drainage."))
+ ((bool, isImbibitionActivated, false,, "Activates imbibition."))
+
,/*TwoPhaseFlowEngineT()*/,
,
.def("getCellIsFictious",&TwoPhaseFlowEngine::cellIsFictious,"Check the connection between pore and boundary. If true, pore throat connects the boundary.")
@@ -166,6 +209,13 @@
.def("setCellHasInterface",&TwoPhaseFlowEngine::setCellHasInterface,"change wheter a cell has a NW-W interface")
.def("savePoreNetwork",&TwoPhaseFlowEngine::savePoreNetwork,"Extract the pore network of the granular material")
.def("getCellLabel",&TwoPhaseFlowEngine::cellLabel,"get cell label. 0 for NW-reservoir; 1 for W-reservoir; others for disconnected W-clusters.")
+ .def("getMinDrainagePc",&TwoPhaseFlowEngine::getMinDrainagePc,"Get the minimum entry capillary pressure for the next drainage step.")
+ .def("getMaxImbibitionPc",&TwoPhaseFlowEngine::getMaxImbibitionPc,"Get the maximum entry capillary pressure for the next imbibition step.")
+ .def("getSaturation",&TwoPhaseFlowEngine::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("invasion",&TwoPhaseFlowEngine::invasion,"Run the drainage invasion.")
+ .def("computeCapillaryForce",&TwoPhaseFlowEngine::computeCapillaryForce,"Compute capillary force. ")
+ .def("saveVtk",&TwoPhaseFlowEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
+ .def("getPotentialPendularSpheresPair",&TwoPhaseFlowEngine::getPotentialPendularSpheresPair,"Get the list of sphere ID pairs of potential pendular liquid bridge.")
)
DECLARE_LOGGER;
=== modified file 'pkg/pfv/UnsaturatedEngine.cpp'
--- pkg/pfv/UnsaturatedEngine.cpp 2016-11-15 18:09:17 +0000
+++ pkg/pfv/UnsaturatedEngine.cpp 2016-11-22 17:29:15 +0000
@@ -23,29 +23,6 @@
double computeCellInterfacialArea(CellHandle cell, int j, double rC);
// void computeSolidLine();
- void computeFacetPoreForcesWithCache(bool onlyCache=false);
- void computeCapillaryForce() {computeFacetPoreForcesWithCache(false);}
-
-
- void invasion();
- ///functions can be shared by two modes
- void invasionSingleCell(CellHandle cell);
- void updatePressure();
- double getMinDrainagePc();
- double getMaxImbibitionPc();
- double getSaturation(bool isSideBoundaryIncluded=false);
- double getSpecificInterfacialArea();
-
- void invasion1();
- void updateReservoirs1();
- void WResRecursion(CellHandle cell);
- void NWResRecursion(CellHandle cell);
- void checkTrap(double pressure);
-
- void invasion2();
- void updateReservoirs2();
-
- void saveVtk(const char* folder) {bool initT=solver->noCache; solver->noCache=false; solver->saveVtk(folder); solver->noCache=initT;}
//record and test functions
void checkLatticeNodeY(double y);
@@ -57,40 +34,19 @@
double getSphericalSubdomainSaturation(Vector3r pos, double radius);
double getCuboidSubdomainSaturation(Vector3r pos1, Vector3r pos2, bool isSideBoundaryIncluded);
double getCuboidSubdomainPorosity(Vector3r pos1, Vector3r pos2, bool isSideBoundaryIncluded);
+ double getSpecificInterfacialArea();
void printSomething();
- boost::python::list getPotentialPendularSpheresPair() {
- RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- boost::python::list bridgeIds;
- FiniteEdgesIterator ed_it = Tri.finite_edges_begin();
- for ( ; ed_it!=Tri.finite_edges_end(); ed_it++ ) {
- if (detectBridge(ed_it)==true) {
- const VertexInfo& vi1=(ed_it->first)->vertex(ed_it->second)->info();
- const VertexInfo& vi2=(ed_it->first)->vertex(ed_it->third)->info();
- const int& id1 = vi1.id();
- const int& id2 = vi2.id();
- bridgeIds.append(boost::python::make_tuple(id1,id2));}}
- return bridgeIds;}
- bool detectBridge(RTriangulation::Finite_edges_iterator& edge);
virtual ~UnsaturatedEngine();
virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(UnsaturatedEngine,TwoPhaseFlowEngine,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(UnsaturatedEngine,TwoPhaseFlowEngine,"Preliminary version engine of a drainage model for unsaturated soils. Note:Air reservoir is on the top; water reservoir is on the bottom.(deprecated engine, use TwoPhaseFlowEngine instead)",
- ((bool, computeForceActivated, true,,"Activate capillary force computation. WARNING: turning off means capillary force is not computed at all, but the drainage can still work."))
((int, windowsNo, 10,, "Number of genrated windows(or zoomed samples)."))
- ((bool, isDrainageActivated, true,, "Activates drainage."))
- ((bool, isImbibitionActivated, false,, "Activates imbibition."))
,,,
- .def("saveVtk",&UnsaturatedEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .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("invasion",&UnsaturatedEngine::invasion,"Run the drainage invasion.")
- .def("computeCapillaryForce",&UnsaturatedEngine::computeCapillaryForce,"Compute capillary force. ")
.def("checkLatticeNodeY",&UnsaturatedEngine::checkLatticeNodeY,(boost::python::arg("y")),"Check the slice of lattice nodes for yNormal(y). 0: out of sphere; 1: inside of sphere.")
.def("getInvadeDepth",&UnsaturatedEngine::getInvadeDepth,"Get NW-phase invasion depth. (the distance from NW-reservoir to front of NW-W interface.)")
.def("getSphericalSubdomainSaturation",&UnsaturatedEngine::getSphericalSubdomainSaturation,(boost::python::arg("pos"),boost::python::arg("radius")),"Get saturation of spherical subdomain defined by (pos, radius). The subdomain exclude boundary pores.")
@@ -100,7 +56,6 @@
.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 temporary function for comparison with experiments, will delete soon")
.def("printSomething",&UnsaturatedEngine::printSomething,"print debug.")
- .def("getPotentialPendularSpheresPair",&UnsaturatedEngine::getPotentialPendularSpheresPair,"Get the list of sphere ID pairs of potential pendular liquid bridge.")
)
DECLARE_LOGGER;
};
@@ -161,310 +116,9 @@
scene->forces.addForce ( V_it->info().id(), force); }}
*/}
-void UnsaturatedEngine::updatePressure()
-{
- boundaryConditions(*solver);
- solver->pressureChanged=true;
- solver->reApplyBoundaryConditions();
- 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) {cell->info().p()=bndCondValue[2];}
- if (cell->info().isNWRes==true) {cell->info().p()=bndCondValue[3];}
- if (isPhaseTrapped) {
- if ( cell->info().isTrapW ) {cell->info().p()=bndCondValue[3]-cell->info().trapCapP;}
- if ( cell->info().isTrapNW) {cell->info().p()=bndCondValue[2]+cell->info().trapCapP;}
- //check cell reservoir info.
- if ( !cell->info().isWRes && !cell->info().isNWRes && !cell->info().isTrapW && !cell->info().isTrapNW ) {cerr<<"ERROR! NOT FIND Cell Info!";}
-// {cell->info().p()=bndCondValue[2]; if (isInvadeBoundary) cerr<<"Something wrong in updatePressure.(isInvadeBoundary)";}
- }
- }
-}
-
-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 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 (cell->info().poreThroatRadius[facet]<0) 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) ) {
- nCell->info().p() = localPressure;
- nCell->info().saturation=localSaturation;
- nCell->info().hasInterface=false;
- if(solver->debugOut) {cerr<<"drainage"<<endl;}
- if (recursiveInvasion) invasionSingleCell(nCell);
- }
-////FIXME:Introduce cell.hasInterface
-// else if( (localPressure-nCell->info().p()>nPcThroat) && (localPressure-nCell->info().p()<nPcBody) && (cell->info().hasInterface==false) && (nCell->info().hasInterface==false) ) {
-// if(solver->debugOut) {cerr<<"invasion paused into pore interface "<<endl;}
-// nCell->info().hasInterface=true;
-// }
-// else continue;
- }
- else if ( (nCell->info().saturation<localSaturation) ) {
- double nPcThroat=surfaceTension/cell->info().poreThroatRadius[facet];
- double nPcBody=surfaceTension/nCell->info().poreBodyRadius;
- if( (nCell->info().p()-localPressure<nPcBody) && (nCell->info().p()-localPressure<nPcThroat) ) {
- nCell->info().p() = localPressure;
- nCell->info().saturation=localSaturation;
- if(solver->debugOut) {cerr<<"imbibition"<<endl;}
- if (recursiveInvasion) invasionSingleCell(nCell);
- }
-//// FIXME:Introduce cell.hasInterface
-// else if ( (nCell->info().p()-localPressure<nPcBody) && (nCell->info().p()-localPressure>nPcThroat) /*&& (cell->info().hasInterface==false) && (nCell->info().hasInterface==false)*/ ) {
-// nCell->info().p() = localPressure;
-// nCell->info().saturation=localSaturation;
-// if(solver->debugOut) {cerr<<"imbibition paused pore interface"<<endl;}
-// nCell->info().hasInterface=true;
-// }
-// else continue;
- }
- else continue;
- }
-}
-///invasion mode 1: withTrap
-void UnsaturatedEngine::invasion1()
-{
- if(solver->debugOut) {cout<<"----start invasion1----"<<endl;}
-
- ///update Pw, Pn according to reservoirInfo.
- updatePressure();
- if(solver->debugOut) {cout<<"----invasion1.updatePressure----"<<endl;}
-
- ///invasionSingleCell by Pressure difference, change Pressure and Saturation.
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- FiniteCellsIterator cellEnd = tri.finite_cells_end();
- if(isDrainageActivated) {
- for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if(cell->info().isNWRes)
- invasionSingleCell(cell);
- }
- }
- if(isImbibitionActivated) {
- for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if(cell->info().isWRes)
- invasionSingleCell(cell);
- }
- }
- if(solver->debugOut) {cout<<"----invasion1.invasionSingleCell----"<<endl;}
-
- ///update W, NW reservoirInfo according to cell->info().saturation
- updateReservoirs1();
- if(solver->debugOut) {cout<<"----invasion1.update W, NW reservoirInfo----"<<endl;}
-
- ///search new trapped W-phase/NW-phase, assign trapCapP, isTrapW/isTrapNW flag for new trapped phases. But at this moment, the new trapped W/NW cells.P= W/NW-Res.P. They will be updated in next updatePressure() func.
- checkTrap(bndCondValue[3]-bndCondValue[2]);
- if(solver->debugOut) {cout<<"----invasion1.checkWTrap----"<<endl;}
-
- ///update trapped W-phase/NW-phase Pressure //FIXME: is this necessary?
- 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;}
-
- if(isCellLabelActivated) updateCellLabel();
- if(solver->debugOut) {cout<<"----update cell labels----"<<endl;}
-}
-
-///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().isFictious) && (!cell->info().Pcondition) && (!isInvadeBoundary) ) continue;
- if( (cell->info().isWRes) || (cell->info().isNWRes) || (cell->info().isTrapW) || (cell->info().isTrapNW) ) continue;
- cell->info().trapCapP=pressure;
- if(cell->info().saturation==1.0) cell->info().isTrapW=true;
- if(cell->info().saturation==0.0) cell->info().isTrapNW=true;
- }
-}
-
-void UnsaturatedEngine::updateReservoirs1()
-{
- 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;
- }
-
- for (FlowSolver::VCellIterator it = solver->boundingCells[2].begin(); it != solver->boundingCells[2].end(); it++) {
- if ((*it)==NULL) continue;
- WResRecursion(*it);
- }
-
- for (FlowSolver::VCellIterator it = solver->boundingCells[3].begin(); it != solver->boundingCells[3].end(); it++) {
- if ((*it)==NULL) continue;
- NWResRecursion(*it);
- }
-}
-
-void UnsaturatedEngine::WResRecursion(CellHandle cell)
-{
- 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;
-// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
- if (nCell->info().saturation != 1.0) continue;
- 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);
- }
-}
-
-void UnsaturatedEngine::NWResRecursion(CellHandle cell)
-{
- 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;
-// if ( (nCell->info().isFictious) && (!isInvadeBoundary) ) continue;
- if (nCell->info().saturation != 0.0) continue;
- 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);
- }
-}
-
-///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, change Pressure and Saturation.
- RTriangulation& tri = solver->T[solver->currentTes].Triangulation();
- FiniteCellsIterator cellEnd = tri.finite_cells_end();
- if(isDrainageActivated) {
- for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if(cell->info().isNWRes)
- invasionSingleCell(cell);
- }
- }
- ///drainageSingleCell by Pressure difference, change Pressure and Saturation.
- if(isImbibitionActivated) {
- for ( FiniteCellsIterator cell = tri.finite_cells_begin(); cell != cellEnd; cell++ ) {
- if(cell->info().isWRes)
- invasionSingleCell(cell);
- }
- }
-
- if(solver->debugOut) {cout<<"----invasion2.invasionSingleCell----"<<endl;}
-
- ///update W, NW reservoirInfo according to Pressure
- updateReservoirs2();
- if(solver->debugOut) {cout<<"----drainage2.update W, NW reservoirInfo----"<<endl;}
-
-}
-
-void UnsaturatedEngine::updateReservoirs2()
-{
- 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[2]) {cell->info().isWRes=true; cell->info().isNWRes=false;}
- else if (cell->info().p()==bndCondValue[3]) {cell->info().isNWRes=true; cell->info().isWRes=false;}
- else {cerr<<"drainage mode2: updateReservoir Error!"<<endl;}
- }
-}
-
-double UnsaturatedEngine::getMinDrainagePc()
-{
- 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().isNWRes == 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().isWRes == true && cell->info().poreThroatRadius[facet]>0) {
- 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) {
- cout << "End drainage !" << endl;
- return nextEntry=0;
- }
- 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 && cell->info().poreThroatRadius[facet]>0) {
- double nCellP = std::min( (surfaceTension/nCell->info().poreBodyRadius), (surfaceTension/cell->info().poreThroatRadius[facet]));
- 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;
- 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()
{
@@ -588,20 +242,6 @@
cerr<<id1<<" "<<id2<<endl;}
}
-bool UnsaturatedEngine::detectBridge(RTriangulation::Finite_edges_iterator& edge)
-{
- bool dryBridgeExist=true;
- const RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- RTriangulation::Cell_circulator cell1 = Tri.incident_cells(*edge);
- RTriangulation::Cell_circulator cell0 = cell1++;
- if(cell0->info().saturation==1) {dryBridgeExist=false; return dryBridgeExist;}
- else {
- while (cell1!=cell0) {
- if(cell1->info().saturation==1) {dryBridgeExist=false;break;}
- else cell1++;}
- return dryBridgeExist;
- }
-}
//----------temporary functions for comparison with experiment-----------------------
void UnsaturatedEngine::initializeCellWindowsID()
@@ -716,156 +356,6 @@
//--------------end of comparison with experiment----------------------------
-///compute forces
-void UnsaturatedEngine::computeFacetPoreForcesWithCache(bool onlyCache)
-{
- RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- CVector nullVect(0,0,0);
- //reset forces
- if (!onlyCache) for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) v->info().forces=nullVect;
-// #ifdef parallel_forces
-// if (solver->noCache) {
-// solver->perVertexUnitForce.clear(); solver->perVertexPressure.clear();
-// solver->perVertexUnitForce.resize(solver->T[solver->currentTes].maxId+1);
-// solver->perVertexPressure.resize(solver->T[solver->currentTes].maxId+1);}
-// #endif
-// CellHandle neighbourCell;
-// VertexHandle mirrorVertex;
- CVector tempVect;
- //FIXME : Ema, be carefull with this (noCache), it needs to be turned true after retriangulation
- if (solver->noCache) {//WARNING:all currentTes must be solver->T[solver->currentTes], should NOT be solver->T[currentTes]
- for (FlowSolver::VCellIterator cellIt=solver->T[solver->currentTes].cellHandles.begin(); cellIt!=solver->T[solver->currentTes].cellHandles.end(); cellIt++) {
- CellHandle& cell = *cellIt;
- //reset cache
- for (int k=0; k<4; k++) cell->info().unitForceVectors[k]=nullVect;
-
- for (int j=0; j<4; j++) if (!Tri.is_infinite(cell->neighbor(j))) {
- const CVector& Surfk = cell->info().facetSurfaces[j];
- //FIXME : later compute that fluidSurf only once in hydraulicRadius, for now keep full surface not modified in cell->info for comparison with other forces schemes
- //The ratio void surface / facet surface
- Real area = sqrt(Surfk.squared_length());
- if (area<=0) cerr <<"AREA <= 0!! AREA="<<area<<endl;
- CVector facetNormal = Surfk/area;
- const std::vector<CVector>& crossSections = cell->info().facetSphereCrossSections;
- CVector fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
- /// handle fictious vertex since we can get the projected surface easily here
- if (cell->vertex(j)->info().isFictious) {
- Real projSurf=std::abs(Surfk[solver->boundary(cell->vertex(j)->info().id()).coordinate]);
- tempVect=-projSurf*solver->boundary(cell->vertex(j)->info().id()).normal;
- cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p();
- //define the cached value for later use with cache*p
- cell->info().unitForceVectors[j]=cell->info().unitForceVectors[j]+ tempVect;
- }
- /// Apply weighted forces f_k=sqRad_k/sumSqRad*f
- CVector facetUnitForce = -fluidSurfk*cell->info().solidLine[j][3];
- CVector facetForce = cell->info().p()*facetUnitForce;
-
- for (int y=0; y<3; y++) {
- cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces + facetForce*cell->info().solidLine[j][y];
- //add to cached value
- cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]+facetUnitForce*cell->info().solidLine[j][y];
- //uncomment to get total force / comment to get only pore tension forces
- if (!cell->vertex(facetVertices[j][y])->info().isFictious) {
- cell->vertex(facetVertices[j][y])->info().forces = cell->vertex(facetVertices[j][y])->info().forces -facetNormal*cell->info().p()*crossSections[j][y];
- //add to cached value
- cell->info().unitForceVectors[facetVertices[j][y]]=cell->info().unitForceVectors[facetVertices[j][y]]-facetNormal*crossSections[j][y];
- }
- }
-// #ifdef parallel_forces
-// solver->perVertexUnitForce[cell->vertex(j)->info().id()].push_back(&(cell->info().unitForceVectors[j]));
-// solver->perVertexPressure[cell->vertex(j)->info().id()].push_back(&(cell->info().p()));
-// #endif
- }
- }
- solver->noCache=false;//cache should always be defined after execution of this function
- }
- if (onlyCache) return;
-
-// else {//use cached values when triangulation doesn't change
-// #ifndef parallel_forces
- for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); cell++) {
- for (int yy=0; yy<4; yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();
- }
-
-// #else
-// #pragma omp parallel for num_threads(ompThreads)
-// for (int vn=0; vn<= solver->T[solver->currentTes].maxId; vn++) {
-// VertexHandle& v = solver->T[solver->currentTes].vertexHandles[vn];
-// const int& id = v->info().id();
-// CVector tf (0,0,0);
-// int k=0;
-// for (vector<const Real*>::iterator c = solver->perVertexPressure[id].begin(); c != solver->perVertexPressure[id].end(); c++)
-// tf = tf + (*(solver->perVertexUnitForce[id][k++]))*(**c);
-// v->info().forces = tf;
-// }
-// #endif
-// }
- if (solver->debugOut) {
- CVector totalForce = nullVect;
- for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
- if (!v->info().isFictious) totalForce = totalForce + v->info().forces;
- else if (solver->boundary(v->info().id()).flowCondition==1) totalForce = totalForce + v->info().forces;
- }
- cout << "totalForce = "<< totalForce << endl;
- }
-}
-
-//#########################################################
-// CONVECTIVE DRYING EXTENSION
-//#########################################################
-
-class PhaseCluster : public Serializable
-{
- double totalCellVolume;
- public :
-
-
- virtual ~PhaseCluster();
- vector<TwoPhaseFlowEngine::CellHandle> pores;
- TwoPhaseFlowEngine::RTriangulation* tri;
-
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PhaseCluster,Serializable,"Preliminary.",
- ((int,label,-1,,"Unique label of this cluster, should be reflected in pores of this cluster."))
- ((double,volume,0,,"cumulated volume of all pores."))
- ((double,entryPc,0,,"smallest entry capillary pressure."))
- ((int,entryPore,0,,"the pore of the cluster incident to the throat with smallest entry Pc."))
- ((double,interfacialArea,0,,"interfacial area of the cluster"))
- ,,,
- )
-};
-
-REGISTER_SERIALIZABLE(PhaseCluster);
-YADE_PLUGIN((PhaseCluster));
-
-PhaseCluster::~PhaseCluster(){}
-
-
-class DryingEngine : public UnsaturatedEngine
-{
- public :
- virtual ~DryingEngine();
- vector<shared_ptr<PhaseCluster> > clusters;
-
- boost::python::list pyClusters() {
- boost::python::list ret;
- for(vector<shared_ptr<PhaseCluster> >::iterator it=clusters.begin(); it!=clusters.end(); ++it) ret.append(*it);
- return ret;}
-
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DryingEngine,UnsaturatedEngine,"Extended TwoPhaseFlowEngine for application to convective drying.",
-// ((shared_ptr<PhaseCluster> , cluster,new PhaseCluster,,"The list of clusters"))
- ,,,
- .def("getClusters",&DryingEngine::pyClusters/*,(boost::python::arg("folder")="./VTK")*/,"Save pressure field in vtk format. Specify a folder name for output.")
- )
- DECLARE_LOGGER;
-};
-
-DryingEngine::~DryingEngine(){};
-
-REGISTER_SERIALIZABLE(DryingEngine);
-YADE_PLUGIN((DryingEngine));
-
-
#endif //TWOPHASEFLOW
#endif //FLOW_ENGINE
