From e356880901512242ae34fa0367fd047f28a2ac90 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20Sch=C3=B6nherr?= <schoen@irmb.tu-bs.de>
Date: Mon, 4 Dec 2017 08:56:42 +0100
Subject: [PATCH] porous media v2 there is one more bug

---
 .../Calculation/PorousMedia.cpp               |  20 +-
 .../Calculation/PorousMedia.h                 |  13 +-
 src/VirtualFluids_GPU/Init/SetParameter.cpp   |   1 +
 src/VirtualFluids_GPU/LBM/LB.h                |   1 +
 src/VirtualFluids_GPU/LBM/Simulation.cpp      | 252 +++++++++++-------
 src/VirtualFluids_GPU/LBM/Simulation.h        |   1 +
 .../Output/UnstructuredGridWriter.hpp         | 143 ++++++++++
 src/VirtualFluids_GPU/Output/WriteData.cpp    |   8 +
 src/VirtualFluids_GPU/Parameter/Parameter.cpp |  30 ++-
 src/VirtualFluids_GPU/Parameter/Parameter.h   |   8 +-
 10 files changed, 366 insertions(+), 111 deletions(-)

diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
index 83e7eb0ae..6c1dd5df6 100644
--- a/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
@@ -19,7 +19,7 @@ PorousMedia::PorousMedia()
 	setCoordinatesToZero();
 }
 
-PorousMedia::PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM)
+PorousMedia::PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM, unsigned int level)
 {
 	this->porosity = porosity;
 	this->geoID = geoID;
@@ -27,8 +27,9 @@ PorousMedia::PorousMedia(double porosity, unsigned int geoID, double darcySI, do
 	this->forchheimerSI = forchheimerSI;
 	this->dxLBM = dxLBM;
 	this->dtLBM = dtLBM;
-	this->forchheimerLBM = this->forchheimerSI * this->dtLBM;
-	this->darcyLBM = this->darcySI * this->dxLBM * this->dxLBM;
+	this->forchheimerLBM = this->forchheimerSI * this->dxLBM;
+	this->darcyLBM = this->darcySI * this->dtLBM;
+	this->level = level;
 	setCoordinatesToZero();
 }
 
@@ -66,6 +67,11 @@ void PorousMedia::setDeviceNodeIDsPM(unsigned int* deviceNodeIDsPM)
 	this->deviceNodeIDsPM = deviceNodeIDsPM;
 }
 
+void PorousMedia::setSizePM(unsigned int sizePM)
+{
+	this->sizePM = sizePM;
+}
+
 //void PorousMedia::definePMarea(Parameter* para, unsigned int level)
 //{
 //	unsigned int counter = 0;
@@ -99,5 +105,13 @@ double PorousMedia::getDarcyLBM() {					return this->darcyLBM; }
 double PorousMedia::getForchheimerLBM() {			return this->forchheimerLBM; }
 unsigned int PorousMedia::getGeoID() {				return this->geoID; }
 unsigned int PorousMedia::getSizePM() {				return this->sizePM; }
+unsigned int PorousMedia::getLevelPM() {			return this->level; }
 unsigned int* PorousMedia::getHostNodeIDsPM() {		return this->hostNodeIDsPM; }
 unsigned int* PorousMedia::getDeviceNodeIDsPM() {	return this->deviceNodeIDsPM; }
+double PorousMedia::getStartX(){					return this->startCoordX; }
+double PorousMedia::getStartY(){					return this->startCoordY; }
+double PorousMedia::getStartZ(){					return this->startCoordZ; }
+double PorousMedia::getEndX(){						return this->endCoordX; }
+double PorousMedia::getEndY(){						return this->endCoordY; }
+double PorousMedia::getEndZ(){						return this->endCoordZ; }
+
diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.h b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
index da0c9741b..a6d49037d 100644
--- a/src/VirtualFluids_GPU/Calculation/PorousMedia.h
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
@@ -14,7 +14,7 @@ class PorousMedia
 {
 public:
 	PorousMedia();
-	PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM);
+	PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM, unsigned int level);
 	~PorousMedia();
 
 	//setter
@@ -26,6 +26,7 @@ public:
 	//void definePMarea(Parameter* para, unsigned int level);
 	void setHostNodeIDsPM(unsigned int* hostNodeIDsPM);
 	void setDeviceNodeIDsPM(unsigned int* deviceNodeIDsPM);
+	void setSizePM(unsigned int sizePM);
 
 	//getter
 	double getPorosity();
@@ -35,8 +36,15 @@ public:
 	double getForchheimerLBM();
 	unsigned int getGeoID();
 	unsigned int getSizePM();
+	unsigned int getLevelPM();
 	unsigned int* getHostNodeIDsPM();
 	unsigned int* getDeviceNodeIDsPM();
+	double getStartX();
+	double getStartY();
+	double getStartZ();
+	double getEndX();
+	double getEndY();
+	double getEndZ();
 
 private:
 	double porosity;
@@ -53,8 +61,9 @@ private:
 	double endCoordY;
 	double endCoordZ;
 	unsigned int geoID;
-	std::vector< unsigned int > nodeIDsPorousMedia;
+	//std::vector< unsigned int > nodeIDsPorousMedia;
 	unsigned int sizePM;
+	unsigned int level;
 	unsigned int *hostNodeIDsPM;
 	unsigned int *deviceNodeIDsPM;
 
diff --git a/src/VirtualFluids_GPU/Init/SetParameter.cpp b/src/VirtualFluids_GPU/Init/SetParameter.cpp
index dc6e15c91..98ac416ca 100644
--- a/src/VirtualFluids_GPU/Init/SetParameter.cpp
+++ b/src/VirtualFluids_GPU/Init/SetParameter.cpp
@@ -35,6 +35,7 @@ void setParameters(Parameter* para, Communicator* comm, std::string &cstr)
    para->setConcFile(            StringUtil::toBool( cf.getValue( "UseConcFile" )));                       
    para->setUseMeasurePoints(    StringUtil::toBool( cf.getValue( "UseMeasurePoints")));
    para->setUseWale(             StringUtil::toBool( cf.getValue( "UseWale" )));
+   para->setSimulatePorousMedia( StringUtil::toBool( cf.getValue( "SimulatePorousMedia" )));
    para->setD3Qxx(               StringUtil::toInt(  cf.getValue( "D3Qxx" )));
    para->setMaxLevel(            StringUtil::toInt(  cf.getValue( "NOGL" )));                             
    para->setTEnd(                StringUtil::toInt(  cf.getValue( "TimeEnd" )));                          
diff --git a/src/VirtualFluids_GPU/LBM/LB.h b/src/VirtualFluids_GPU/LBM/LB.h
index e454f8468..ed99d1c74 100644
--- a/src/VirtualFluids_GPU/LBM/LB.h
+++ b/src/VirtualFluids_GPU/LBM/LB.h
@@ -102,6 +102,7 @@ typedef struct InitCond{
    doubflo delta_rho, delta_press;
    bool  printFiles, readGeo, doRestart, doCheckPoint, isGeo, isProp, isCp, calcMedian, GeometryValues, isConc, is2ndOrderMoments, is3rdOrderMoments, isHighOrderMoments, isWale, isMeasurePoints;
    bool isGeoNormal, isInflowNormal, isOutflowNormal;
+   bool simulatePorousMedia;
 } InitCondition;
 
 //Interface Cells
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.cpp b/src/VirtualFluids_GPU/LBM/Simulation.cpp
index 1569eda25..debcfc436 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -183,6 +183,16 @@ void Simulation::init(std::string &cstr)
 	   output << "done.\n";
    }
 
+   //////////////////////////////////////////////////////////////////////////
+   //Porous Media
+   //////////////////////////////////////////////////////////////////////////
+   if (para->getSimulatePorousMedia())
+   {
+	   output << "define area(s) of porous media...";
+	   porousMedia();
+	   output << "done.\n";
+   }
+
    //////////////////////////////////////////////////////////////////////////
    //enSightGold
    //////////////////////////////////////////////////////////////////////////
@@ -545,100 +555,105 @@ void Simulation::run()
 		 //getLastCudaError("KernelCasSPKum27 execution failed");
 		 //////////////////////////////////////////////////////////////////////////
 
-		////////////////////////////////////////////////////////////////////////////
-		////Wale 
-		//if (para->getUseWale())
-		//{
-		//	KernelWaleCumOneCompSP27(para->getParD(0)->numberofthreads,
-		//							 para->getParD(0)->omega,			
-		//							 para->getParD(0)->geoSP, 
-		//							 para->getParD(0)->neighborX_SP, 
-		//							 para->getParD(0)->neighborY_SP, 
-		//							 para->getParD(0)->neighborZ_SP,
-		//							 para->getParD(0)->neighborWSB_SP,
-		//						     para->getParD(0)->vx_SP,        
-		//						     para->getParD(0)->vy_SP,        
-		//						     para->getParD(0)->vz_SP,        
-		//							 para->getParD(0)->d0SP.f[0],
-		//							 para->getParD(0)->turbViscosity,
-		//							 para->getParD(0)->size_Mat_SP,
-		//							 para->getParD(0)->size_Array_SP,
-		//							 0,
-		//							 para->getForcesDev(),
-		//							 para->getParD(0)->evenOrOdd); 
-		//	getLastCudaError("KernelWaleCumOneCompSP27 execution failed");
-
-		//	//KernelWaleCumAA2016CompSP27(para->getParD(0)->numberofthreads,
-		//	//							para->getParD(0)->omega,			
-		//	//							para->getParD(0)->geoSP, 
-		//	//							para->getParD(0)->neighborX_SP, 
-		//	//							para->getParD(0)->neighborY_SP, 
-		//	//							para->getParD(0)->neighborZ_SP,
-		//	//							para->getParD(0)->neighborWSB_SP,
-		//	//							para->getParD(0)->vx_SP,        
-		//	//							para->getParD(0)->vy_SP,        
-		//	//							para->getParD(0)->vz_SP,        
-		//	//							para->getParD(0)->d0SP.f[0],
-		//	//							para->getParD(0)->turbViscosity,
-		//	//							para->getParD(0)->size_Mat_SP,
-		//	//							para->getParD(0)->size_Array_SP,
-		//	//							0,
-		//	//							para->getForcesDev(),
-		//	//							para->getParD(0)->evenOrOdd); 
-		//	//getLastCudaError("KernelWaleCumAA2016CompSP27 execution failed");
-
-		//} 
-		//else
-		//{
-		//	KernelKumNewCompSP27(para->getParD(0)->numberofthreads,       
-		//						 para->getParD(0)->omega,			
-		//						 para->getParD(0)->geoSP, 
-		//						 para->getParD(0)->neighborX_SP, 
-		//						 para->getParD(0)->neighborY_SP, 
-		//						 para->getParD(0)->neighborZ_SP,
-		//						 para->getParD(0)->d0SP.f[0],    
-		//						 para->getParD(0)->size_Mat_SP,
-		//						 para->getParD(0)->size_Array_SP,
-		//						 0,
-		//						 para->getForcesDev(),
-		//						 para->getParD(0)->evenOrOdd); 
-		//	getLastCudaError("KernelCasSPKum27 execution failed");
-
-		//	//KernelKumAA2016CompSP27(para->getParD(0)->numberofthreads,       
-		//	//					 para->getParD(0)->omega, 
-		//	//					 para->getParD(0)->geoSP, 
-		//	//					 para->getParD(0)->neighborX_SP, 
-		//	//					 para->getParD(0)->neighborY_SP, 
-		//	//					 para->getParD(0)->neighborZ_SP,
-		//	//					 para->getParD(0)->d0SP.f[0],    
-		//	//					 para->getParD(0)->size_Mat_SP,
-		//	//					 0,
-		//	//					 para->getForcesDev(),
-		//	//					 para->getParD(0)->evenOrOdd); 
-		//	//getLastCudaError("KernelKumAA2016CompSP27 execution failed");
-
-		//}
-		////////////////////////////////////////////////////////////////////////////
+		//////////////////////////////////////////////////////////////////////////
+		//Wale 
+		if (para->getUseWale())
+		{
+			KernelWaleCumOneCompSP27(para->getParD(0)->numberofthreads,
+									 para->getParD(0)->omega,			
+									 para->getParD(0)->geoSP, 
+									 para->getParD(0)->neighborX_SP, 
+									 para->getParD(0)->neighborY_SP, 
+									 para->getParD(0)->neighborZ_SP,
+									 para->getParD(0)->neighborWSB_SP,
+								     para->getParD(0)->vx_SP,        
+								     para->getParD(0)->vy_SP,        
+								     para->getParD(0)->vz_SP,        
+									 para->getParD(0)->d0SP.f[0],
+									 para->getParD(0)->turbViscosity,
+									 para->getParD(0)->size_Mat_SP,
+									 para->getParD(0)->size_Array_SP,
+									 0,
+									 para->getForcesDev(),
+									 para->getParD(0)->evenOrOdd); 
+			getLastCudaError("KernelWaleCumOneCompSP27 execution failed");
+
+			//KernelWaleCumAA2016CompSP27(para->getParD(0)->numberofthreads,
+			//							para->getParD(0)->omega,			
+			//							para->getParD(0)->geoSP, 
+			//							para->getParD(0)->neighborX_SP, 
+			//							para->getParD(0)->neighborY_SP, 
+			//							para->getParD(0)->neighborZ_SP,
+			//							para->getParD(0)->neighborWSB_SP,
+			//							para->getParD(0)->vx_SP,        
+			//							para->getParD(0)->vy_SP,        
+			//							para->getParD(0)->vz_SP,        
+			//							para->getParD(0)->d0SP.f[0],
+			//							para->getParD(0)->turbViscosity,
+			//							para->getParD(0)->size_Mat_SP,
+			//							para->getParD(0)->size_Array_SP,
+			//							0,
+			//							para->getForcesDev(),
+			//							para->getParD(0)->evenOrOdd); 
+			//getLastCudaError("KernelWaleCumAA2016CompSP27 execution failed");
+
+		} 
+		else
+		{
+			KernelKumNewCompSP27(para->getParD(0)->numberofthreads,       
+								 para->getParD(0)->omega,			
+								 para->getParD(0)->geoSP, 
+								 para->getParD(0)->neighborX_SP, 
+								 para->getParD(0)->neighborY_SP, 
+								 para->getParD(0)->neighborZ_SP,
+								 para->getParD(0)->d0SP.f[0],    
+								 para->getParD(0)->size_Mat_SP,
+								 para->getParD(0)->size_Array_SP,
+								 0,
+								 para->getForcesDev(),
+								 para->getParD(0)->evenOrOdd); 
+			getLastCudaError("KernelCasSPKum27 execution failed");
+
+			//KernelKumAA2016CompSP27(para->getParD(0)->numberofthreads,       
+			//					 para->getParD(0)->omega, 
+			//					 para->getParD(0)->geoSP, 
+			//					 para->getParD(0)->neighborX_SP, 
+			//					 para->getParD(0)->neighborY_SP, 
+			//					 para->getParD(0)->neighborZ_SP,
+			//					 para->getParD(0)->d0SP.f[0],    
+			//					 para->getParD(0)->size_Mat_SP,
+			//					 0,
+			//					 para->getForcesDev(),
+			//					 para->getParD(0)->evenOrOdd); 
+			//getLastCudaError("KernelKumAA2016CompSP27 execution failed");
+
+		}
+		//////////////////////////////////////////////////////////////////////////
 
 
 		//////////////////////////////////////////////////////////////////////////
 		//porous media
-		KernelPMCumOneCompSP27( para->getParD(0)->numberofthreads,
-								para->getParD(0)->omega,			
-								para->getParD(0)->geoSP, 
-								para->getParD(0)->neighborX_SP, 
-								para->getParD(0)->neighborY_SP, 
-								para->getParD(0)->neighborZ_SP,
-								para->getParD(0)->d0SP.f[0],    
-								para->getParD(0)->size_Mat_SP,
-								0,
-								para->getForcesDev(),
-								0.5,
-								134.4,
-								0.0068287,
-								GEO_PM_1,
-								para->getParD(0)->evenOrOdd); 
-		getLastCudaError("KernelPMCumOneCompSP27 execution failed");
+		if (para->getSimulatePorousMedia())
+		{
+			KernelPMCumOneCompSP27( para->getParD(0)->numberofthreads,
+									para->getParD(0)->omega,			
+									para->getParD(0)->geoSP, 
+									para->getParD(0)->neighborX_SP, 
+									para->getParD(0)->neighborY_SP, 
+									para->getParD(0)->neighborZ_SP,
+									para->getParD(0)->d0SP.f[0],    
+									para->getParD(0)->size_Mat_SP,
+									0,
+									para->getForcesDev(),
+									pm0->getPorosity(),
+									pm0->getDarcyLBM(),
+									pm0->getForchheimerLBM(),
+									pm0->getGeoID(),
+									para->getParD(0)->evenOrOdd); 
+			getLastCudaError("KernelPMCumOneCompSP27 execution failed");
+		}
+		//////////////////////////////////////////////////////////////////////////
+
 
 
 
@@ -2608,7 +2623,60 @@ void Simulation::run()
 void Simulation::porousMedia()
 {
 	//Kondensator = porous media 0
-	double porosity = 0.5;
-	//pm0 = new PorousMedia(porosity, GEO_PM_0,)
+	double porosity = 0.7;
+	double darcySI = 137.36; //[1/s]
+	double forchheimerSI = 1037.8; //[1/m]
+	double dxLBM = 0.0004;
+	double dtLBM = 0.000007;
+	unsigned int level = 0;
+	output << "\nnew instance of PM \n";
+	pm0 = new PorousMedia(porosity, GEO_PM_0, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+	output << "setStartCoordinates \n";
+	pm0->setStartCoordinates(20, 0, 0);
+	output << "setEndCoordinates \n";
+	pm0->setEndCoordinates(40, 20, 20);
+	output << "definePMarea \n";
+	definePMarea(pm0);
+}
+void Simulation::definePMarea(PorousMedia* pm)
+{
+	unsigned int counter = 0;
+	unsigned int level = pm->getLevelPM();
+	std::vector< unsigned int > nodeIDsPorousMedia;
+	output << "definePMarea....find nodes \n";
+
+	for (unsigned int i = 0; i < para->getParH(level)->size_Mat_SP; i++)
+	{
+		//output << "definePMarea....find nodes...for loop \n";
+		if (((para->getParH(level)->coordX_SP[i] >= pm->getStartX()) && (para->getParH(level)->coordX_SP[i] <= pm->getEndX())) &&
+			((para->getParH(level)->coordY_SP[i] >= pm->getStartY()) && (para->getParH(level)->coordY_SP[i] <= pm->getEndY())) &&
+			((para->getParH(level)->coordZ_SP[i] >= pm->getStartZ()) && (para->getParH(level)->coordZ_SP[i] <= pm->getEndZ())) )
+		{
+			//output << "definePMarea....find nodes...if area \n";
+			if (para->getParH(level)->geoSP[i] >= GEO_FLUID)
+			{
+				//output << "definePMarea....find nodes...if area and fluid \n";
+				para->getParH(level)->geoSP[i] = pm->getGeoID();
+				nodeIDsPorousMedia.push_back(i);
+				counter++;
+			}
+		}
+	}
+
+	output << "definePMarea....cuda alloc \n";
+	pm->setSizePM(counter);
+	para->cudaAllocPorousMedia(pm0, level);
+	unsigned int *tpmArrayIDs = pm->getHostNodeIDsPM();
+
+	output << "definePMarea....copy vector to array \n";
+	for (unsigned int j = 0; j <= pm->getSizePM(); j++)
+	{
+		tpmArrayIDs[j] = nodeIDsPorousMedia[j];
+	}
+
+	output << "definePMarea....setHostNodeIDsPM \n";
+	pm->setHostNodeIDsPM(tpmArrayIDs);
 
+	output << "definePMarea....cudaCopyPorousMedia \n";
+	para->cudaCopyPorousMedia(pm, level);
 }
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.h b/src/VirtualFluids_GPU/LBM/Simulation.h
index 7fc3539d9..ec77b6c9e 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.h
+++ b/src/VirtualFluids_GPU/LBM/Simulation.h
@@ -24,6 +24,7 @@ public:
 	void init(std::string &cstr);
 	void bulk();
 	void porousMedia();
+	void definePMarea(PorousMedia* pm);
 protected:
 
 	Buffer2D <doubflo> sbuf_t; 
diff --git a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
index 4284a1792..bc3089bf6 100644
--- a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
+++ b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
@@ -442,6 +442,149 @@ namespace UnstrucuredGridWriter
 
 
 
+	//////////////////////////////////////////////////////////////////////////
+	void writeUnstrucuredGridPM(Parameter* para, int level, vector<string >& fname)
+	{
+		vector< UbTupleFloat3 > nodes;
+		vector< UbTupleUInt8 > cells;
+		//vector< UbTupleUInt8 > cells2;
+		vector< string > nodedatanames;
+		nodedatanames.push_back("press");
+		nodedatanames.push_back("rho");
+		nodedatanames.push_back("vx1");
+		nodedatanames.push_back("vx2");
+		nodedatanames.push_back("vx3");
+		nodedatanames.push_back("geo");
+		unsigned int number1, number2, number3, number4, number5, number6, number7, number8;
+		unsigned int dn1, dn2, dn3, dn4, dn5, dn6, dn7, dn8;
+		bool neighborsFluid;
+		double vxmax = 0;
+		unsigned int startpos = 0;
+		unsigned int endpos = 0;
+		unsigned int sizeOfNodes = 0;
+		vector< vector< double > > nodedata(nodedatanames.size());
+
+
+		//printf("\n test for if... \n");
+		for (unsigned int part = 0; part < fname.size(); part++)
+		{
+			vxmax = 0;
+			//printf("\n test in if I... \n");
+			//////////////////////////////////////////////////////////////////////////
+			if (((part + 1)*para->getlimitOfNodesForVTK()) > para->getParH(level)->size_Mat_SP)
+			{
+				sizeOfNodes = para->getParH(level)->size_Mat_SP - (part * para->getlimitOfNodesForVTK());
+			}
+			else
+			{
+				sizeOfNodes = para->getlimitOfNodesForVTK();
+			}
+			//////////////////////////////////////////////////////////////////////////
+			startpos = part * para->getlimitOfNodesForVTK();
+			endpos = startpos + sizeOfNodes;
+			//////////////////////////////////////////////////////////////////////////
+			cells.clear();
+			nodes.resize(sizeOfNodes);
+			nodedata[0].resize(sizeOfNodes);
+			nodedata[1].resize(sizeOfNodes);
+			nodedata[2].resize(sizeOfNodes);
+			nodedata[3].resize(sizeOfNodes);
+			nodedata[4].resize(sizeOfNodes);
+			nodedata[5].resize(sizeOfNodes);
+			//////////////////////////////////////////////////////////////////////////
+			//int counter = 0;
+			//////////////////////////////////////////////////////////////////////////
+			//printf("\n test in if II... \n");
+
+			for (unsigned int pos = startpos; pos < endpos; pos++)
+			{
+				if (/*para->getParH(level)->geoSP[pos] >= GEO_FLUID*/true)
+				{
+					//////////////////////////////////////////////////////////////////////////
+					double x1 = para->getParH(level)->coordX_SP[pos];
+					double x2 = para->getParH(level)->coordY_SP[pos];
+					double x3 = para->getParH(level)->coordZ_SP[pos];
+					//////////////////////////////////////////////////////////////////////////
+					number1 = pos;
+					dn1 = pos - startpos;
+					neighborsFluid = true;
+					//////////////////////////////////////////////////////////////////////////
+					//printf("\n test vor node data... \n");
+					nodes[dn1] = (makeUbTuple((float)(x1), (float)(x2), (float)(x3)));
+					nodedata[0][dn1] = (double)para->getParH(level)->press_SP[pos] / (double)3.0 * (double)para->getDensityRatio() * (double)para->getVelocityRatio() * (double)para->getVelocityRatio();
+					nodedata[1][dn1] = (double)para->getParH(level)->rho_SP[pos] / (double)3.0 * (double)para->getDensityRatio() * (double)para->getVelocityRatio() * (double)para->getVelocityRatio();
+					nodedata[2][dn1] = (double)para->getParH(level)->vx_SP[pos] * (double)para->getVelocityRatio();
+					nodedata[3][dn1] = (double)para->getParH(level)->vy_SP[pos] * (double)para->getVelocityRatio();
+					nodedata[4][dn1] = (double)para->getParH(level)->vz_SP[pos] * (double)para->getVelocityRatio();
+					nodedata[5][dn1] = (double)para->getParH(level)->geoSP[pos];
+					//////////////////////////////////////////////////////////////////////////
+					//printf("\n test vor numbers... \n");
+					number2 = para->getParH(level)->neighborX_SP[number1];
+					number3 = para->getParH(level)->neighborY_SP[number2];
+					number4 = para->getParH(level)->neighborY_SP[number1];
+					number5 = para->getParH(level)->neighborZ_SP[number1];
+					number6 = para->getParH(level)->neighborZ_SP[number2];
+					number7 = para->getParH(level)->neighborZ_SP[number3];
+					number8 = para->getParH(level)->neighborZ_SP[number4];
+					//////////////////////////////////////////////////////////////////////////
+					//printf("\n test vor neighborsFluid... \n");
+					//if (para->getParH(level)->geoSP[number2] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number3] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number4] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number5] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number6] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number7] < GEO_FLUID ||
+					//	para->getParH(level)->geoSP[number8] < GEO_FLUID)  neighborsFluid = false;
+					//////////////////////////////////////////////////////////////////////////
+					//if(neighborsFluid==false) counter++;
+					//////////////////////////////////////////////////////////////////////////
+					//printf("\n test vor numbers and neighborsFluid... \n");
+					if (number2 > endpos ||
+						number3 > endpos ||
+						number4 > endpos ||
+						number5 > endpos ||
+						number6 > endpos ||
+						number7 > endpos ||
+						number8 > endpos)  neighborsFluid = false;
+					//////////////////////////////////////////////////////////////////////////
+					//if(neighborsFluid==false) counter++;
+					//////////////////////////////////////////////////////////////////////////
+					//printf("\n test vor dn... \n");
+					dn2 = number2 - startpos;
+					dn3 = number3 - startpos;
+					dn4 = number4 - startpos;
+					dn5 = number5 - startpos;
+					dn6 = number6 - startpos;
+					dn7 = number7 - startpos;
+					dn8 = number8 - startpos;
+					//////////////////////////////////////////////////////////////////////////
+					//if( std::fabs(nodedata[2][dn1]) > std::fabs(vxmax) ) vxmax = nodedata[2][dn1];
+					//////////////////////////////////////////////////////////////////////////
+
+					if (isPeriodicCell(para, level, number2, number1, number3, number5))
+						continue;
+
+					//counter++;
+					if (neighborsFluid) cells.push_back(makeUbTuple(dn1, dn2, dn3, dn4, dn5, dn6, dn7, dn8));
+					//////////////////////////////////////////////////////////////////////////
+				}
+				//printf("\n test II... \n");
+			}
+			//printf("\n number of cells: %d at level: %d\n", cells.size(), level);
+			WbWriterVtkXmlBinary::getInstance()->writeOctsWithNodeData(fname[part], nodes, cells, nodedatanames, nodedata);
+			//WbWriterVtkXmlBinary::getInstance()->writeNodesWithNodeData(fname[part], nodes, nodedatanames, nodedata);
+			//////////////////////////////////////////////////////////////////////////
+			//printf("\n vx max: %.1f at level: %d\n", vxmax, level);
+			//printf("\n counter: %d at level: %d\n", counter, level);
+		}
+	}
+	//////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
 	//////////////////////////////////////////////////////////////////////////
 	void writeUnstrucuredGridLTConc(Parameter* para, int level, vector<string >& fname) 
 	{
diff --git a/src/VirtualFluids_GPU/Output/WriteData.cpp b/src/VirtualFluids_GPU/Output/WriteData.cpp
index 715c6f217..37a0ad09b 100644
--- a/src/VirtualFluids_GPU/Output/WriteData.cpp
+++ b/src/VirtualFluids_GPU/Output/WriteData.cpp
@@ -80,6 +80,10 @@ void writeInit(Parameter* para)
 			{
 				UnstrucuredGridWriter::writeUnstrucuredGridLTwithTurbulentViscosity(para, lev, fname);
 			}
+			else if (para->getSimulatePorousMedia())
+			{
+				UnstrucuredGridWriter::writeUnstrucuredGridPM(para, lev, fname);
+			}
 			else
 			{
 				UnstrucuredGridWriter::writeUnstrucuredGridLT(para, lev, fname);
@@ -246,6 +250,10 @@ void writeTimestep(Parameter* para, unsigned int t)
 			{
 				UnstrucuredGridWriter::writeUnstrucuredGridLTwithTurbulentViscosity(para, lev, fname);
 			}
+			else if (para->getSimulatePorousMedia())
+			{
+				UnstrucuredGridWriter::writeUnstrucuredGridPM(para, lev, fname);
+			}
 			else
 			{
 				UnstrucuredGridWriter::writeUnstrucuredGridLT(para, lev, fname);
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
index bca9b54d3..78485cd49 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -1967,9 +1967,9 @@ void Parameter::cudaAllocRandomValues()
 }
 //////////////////////////////////////////////////////////////////////////
 //porous media
-void Parameter::cudaAllocPorousMedia(PorousMedia* pm, int pmID, int lev)
+void Parameter::cudaAllocPorousMedia(PorousMedia* pm, int lev)
 {
-	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm[pmID].getSizePM();
+	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm->getSizePM();
 	unsigned int *tmpIDHost, *tmpIDDevice;
 
 	//Host
@@ -1979,25 +1979,25 @@ void Parameter::cudaAllocPorousMedia(PorousMedia* pm, int pmID, int lev)
 	checkCudaErrors(cudaMalloc((void**) &(tmpIDDevice), mem_size_IDsPM));
 
 	//////////////////////////////////////////////////////////////////////////
-	pm[pmID].setHostNodeIDsPM(tmpIDHost);
-	pm[pmID].setDeviceNodeIDsPM(tmpIDDevice);
+	pm->setHostNodeIDsPM(tmpIDHost);
+	pm->setDeviceNodeIDsPM(tmpIDDevice);
 	//////////////////////////////////////////////////////////////////////////
 	double tmp = (double)mem_size_IDsPM;
 	setMemsizeGPU(tmp, false);
 }
-void Parameter::cudaCopyPorousMedia(PorousMedia* pm, int pmID, int lev)
+void Parameter::cudaCopyPorousMedia(PorousMedia* pm, int lev)
 {
-	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm[pmID].getSizePM();
-	unsigned int *tmpIDHost   = pm[pmID].getHostNodeIDsPM();
-	unsigned int *tmpIDDevice = pm[pmID].getDeviceNodeIDsPM();
+	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm->getSizePM();
+	unsigned int *tmpIDHost   = pm->getHostNodeIDsPM();
+	unsigned int *tmpIDDevice = pm->getDeviceNodeIDsPM();
 	//////////////////////////////////////////////////////////////////////////
 	checkCudaErrors(cudaMemcpy(tmpIDDevice, tmpIDHost, mem_size_IDsPM, cudaMemcpyHostToDevice));
 	//////////////////////////////////////////////////////////////////////////
-	pm[pmID].setDeviceNodeIDsPM(tmpIDDevice);
+	pm->setDeviceNodeIDsPM(tmpIDDevice);
 }
-void Parameter::cudaFreePorousMedia(PorousMedia* pm, int pmID, int lev)
+void Parameter::cudaFreePorousMedia(PorousMedia* pm, int lev)
 {
-	checkCudaErrors(cudaFreeHost(pm[pmID].getHostNodeIDsPM()));
+	checkCudaErrors(cudaFreeHost(pm->getHostNodeIDsPM()));
 }
 //////////////////////////////////////////////////////////////////////////
 //advection diffusion
@@ -2871,6 +2871,10 @@ void Parameter::setUseWale(bool useWale)
 {
 	ic.isWale = useWale;
 }
+void Parameter::setSimulatePorousMedia(bool simulatePorousMedia)
+{
+	ic.simulatePorousMedia = simulatePorousMedia;
+}
 void Parameter::setGridX(std::vector<int> GridX)
 {
 	ic.GridX = GridX;
@@ -4144,6 +4148,10 @@ bool Parameter::getUseWale()
 {
 	return ic.isWale;
 }
+bool Parameter::getSimulatePorousMedia()
+{
+	return ic.simulatePorousMedia;
+}
 bool Parameter::getIsGeometryValues()
 {
 	return ic.GeometryValues;
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.h b/src/VirtualFluids_GPU/Parameter/Parameter.h
index 63c3f2d01..dfaccd56d 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.h
@@ -453,9 +453,9 @@ public:
 
 	//////////////////////////////////////////////////////////////////////////
 	//Porous Media
-	void cudaAllocPorousMedia(PorousMedia* pm, int pmID, int lev);
-	void cudaCopyPorousMedia(PorousMedia* pm, int pmID, int lev);
-	void cudaFreePorousMedia(PorousMedia* pm, int pmID, int lev);
+	void cudaAllocPorousMedia(PorousMedia* pm, int lev);
+	void cudaCopyPorousMedia(PorousMedia* pm, int lev);
+	void cudaFreePorousMedia(PorousMedia* pm, int lev);
 	//////////////////////////////////////////////////////////////////////////
 
 
@@ -662,6 +662,7 @@ public:
 	void setConcFile(bool concFile);
 	void setUseMeasurePoints(bool useMeasurePoints);
 	void setUseWale(bool useWale);
+	void setSimulatePorousMedia(bool simulatePorousMedia);
 	void setclockCycleForMP(doubflo clockCycleForMP);
 	void setDevices(std::vector<int> devices);
 	void setGridX(std::vector<int> GridX);
@@ -890,6 +891,7 @@ public:
 	bool getConcFile();
 	bool getUseMeasurePoints();
 	bool getUseWale();
+	bool getSimulatePorousMedia();
 	double getMemsizeGPU();
 	//1D domain decomposition
 	std::vector<std::string> getPossNeighborFiles(std::string sor);
-- 
GitLab