yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3522
committer: T Sweijen <thomasje100@xxxxxxxxxxx>
timestamp: Tue 2014-11-04 10:40:01 +0100
message:
  update for dynamic two-phase flow
modified:
  pkg/pfv/TwoPhaseFlowEngine.cpp
  pkg/pfv/TwoPhaseFlowEngine.hpp


--
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-11-03 18:59:34 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.cpp	2014-11-04 09:40:01 +0000
@@ -40,21 +40,26 @@
     }
 }
 
-void TwoPhaseFlowEngine:: computePoreSatAtInterface(CellHandle cell)
+double TwoPhaseFlowEngine:: computePoreSatAtInterface(int ID)
 {
     //This function calculates the new saturation of pore at the interface between wetting/nonwetting 
     //filled pores. It substracts the outgoing flux from the water volume  
+    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().id == ID){
+    
     double qout = 0.0, Vw = 0.0;
     
     for(unsigned int ngb = 0; ngb < 4; ngb++)
     {
       //find out/influx of water
-      if(cell->neighbor(ngb)->info().isWRes){
-	qout= qout + std::abs(cell->info().kNorm() [ngb])*(cell->neighbor ( ngb )->info().p()-cell->info().p());    
-      }
+       if(cell->neighbor(ngb)->info().isWRes){
+	qout= qout + std::abs(cell->info().kNorm() [ngb])*(cell->info().p()-cell->neighbor ( ngb )->info().p()); 
+       }
     }
    
-    Vw = cell->info().saturation * cell->info().poreBodyVolume - (qout * scene->dt);  
+    Vw = (cell->info().saturation * cell->info().poreBodyVolume) - (qout * scene->dt);  
     cell->info().saturation = Vw / cell->info().poreBodyVolume;
    
    
@@ -67,6 +72,9 @@
     if(cell->info().saturation > 1.0){
       cout << endl <<"dt was too large!,saturation larger than 1 in cell " << cell->info().id;
       cell->info().saturation = 1.0;}
+      return cell->info().saturation;
+      }  
+}
 }
 
 void TwoPhaseFlowEngine:: computePoreCapillaryPressure(CellHandle cell)
@@ -380,7 +388,7 @@
     return deltaF;
 }
 
-void TwoPhaseFlowEngine::savePhaseVtk(const char* folder)
+void TwoPhaseFlowEngine::savePhaseVtk1(const char* folder)
 {
 // 	RTriangulation& Tri = T[solver->noCache?(!currentTes):currentTes].Triangulation();
 	RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();

=== modified file 'pkg/pfv/TwoPhaseFlowEngine.hpp'
--- pkg/pfv/TwoPhaseFlowEngine.hpp	2014-11-03 18:59:34 +0000
+++ pkg/pfv/TwoPhaseFlowEngine.hpp	2014-11-04 09:40:01 +0000
@@ -13,8 +13,8 @@
 
 //keep this #ifdef as long as you don't really want to realize a final version publicly, it will save compilation time for everyone else
 //when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
-// #define TWOPHASEFLOW
-#ifdef TWOPHASEFLOW
+#define TWOPHASEFLOW
+//#ifdef TWOPHASEFLOW
 
 #include "FlowEngine_TwoPhaseFlowEngineT.hpp"
 
@@ -72,9 +72,9 @@
 	void computePoreBodyVolume();	
 	void computePoreBodyRadius();
 	void computePoreThroatCircleRadius();
-	void computePoreSatAtInterface(CellHandle cell);
+	double computePoreSatAtInterface(int ID);
 	void computePoreCapillaryPressure(CellHandle cell);
-	void savePhaseVtk(const char* folder);
+	void savePhaseVtk1(const char* folder);
 	void computePoreThroatRadius();
 	double computeEffPoreThroatRadius(CellHandle cell, int j);
 	double computeEffPoreThroatRadiusFine(CellHandle cell, int j);
@@ -86,14 +86,33 @@
 	void initializeReservoirs();
 	
 	
+	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;}
+	  for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->info().poreThroatRadius[i]);
+	return ids;
+	}
+
+	boost::python::list getNeighbors(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;}
+	  for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->neighbor(i)->info().id);
+	return ids;
+	}
+	
+	
 	//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
 	void computeOnePhaseFlow() {scene = Omega::instance().getScene().get(); if (!solver) cerr<<"no solver!"<<endl; solver->gaussSeidel(scene->dt);}
 	
 	CELL_SCALAR_GETTER(bool,.isWRes,cellIsWRes)
 	CELL_SCALAR_GETTER(bool,.isNWRes,cellIsNWRes)
+	CELL_SCALAR_SETTER(bool,.isNWRes,setCellIsNWRes)
 	CELL_SCALAR_GETTER(Real,.saturation,cellSaturation)
 	CELL_SCALAR_SETTER(Real,.saturation,setCellSaturation)
+	CELL_SCALAR_GETTER(int,.isFictious,cellIsFictious) //Temporary function to allow for simulations in Python
+	CELL_SCALAR_GETTER(bool,.hasInterface,cellHasInterface) //Temporary function to allow for simulations in Python
+	CELL_SCALAR_SETTER(bool,.hasInterface,setCellHasInterface) //Temporary function to allow for simulations in Python
 
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(TwoPhaseFlowEngine,TwoPhaseFlowEngineT,"documentation here",
 	((double,surfaceTension,0.0728,,"Water Surface Tension in contact with air at 20 Degrees Celsius is: 0.0728(N/m)"))
@@ -104,14 +123,20 @@
 
 	,/*TwoPhaseFlowEngineT()*/,
 	,
-	.def("savePhaseVtk",&TwoPhaseFlowEngine::savePhaseVtk,(boost::python::arg("folder")="./phaseVtk"),"Save the saturation of local pores in vtk format. Sw(NW-pore)=0, Sw(W-pore)=1. Specify a folder name for output.")
+	.def("getCellIsFictious",&TwoPhaseFlowEngine::cellIsFictious,"get an integer to indicate which boundary this is allocated to")
+	.def("setCellIsNWRes",&TwoPhaseFlowEngine::setCellIsNWRes,"set status whether 'wetting reservoir' state")
+	.def("savePhaseVtk1",&TwoPhaseFlowEngine::savePhaseVtk1,(boost::python::arg("folder")="./phaseVtk"),"Save the saturation of local pores in vtk format. Sw(NW-pore)=0, Sw(W-pore)=1. Specify a folder name for output.")
 	.def("getCellIsWRes",&TwoPhaseFlowEngine::cellIsWRes,"get status wrt 'wetting reservoir' state")
 	.def("getCellIsNWRes",&TwoPhaseFlowEngine::cellIsNWRes,"get status wrt 'non-wetting reservoir' state")
 	.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. ")
-
+	.def("computePoreSatAtInterface",&TwoPhaseFlowEngine::computePoreSatAtInterface,(boost::python::arg("ID")),"compute pressure and fluxes in the W-phase")
+	.def("getPoreThroatRadius",&TwoPhaseFlowEngine::cellporeThroatRadius,"get 4 pore throat radii")
+	.def("getNeighbors",&TwoPhaseFlowEngine::getNeighbors,"get 4 neigboring cells")
+	.def("getCellHasInterface",&TwoPhaseFlowEngine::cellHasInterface,"indicates whether a NW-W interface is present within the cell")
+	.def("setCellHasInterface",&TwoPhaseFlowEngine::setCellHasInterface,"change wheter a cell has a NW-W interface")
 	)
 	DECLARE_LOGGER;
 };