Merge pull request #50 from hamidrezanorouzi/main

diffusion model for porosity is added, base class porosity is modified

phasicFlowCoupling/Make/files

@@ -12,6 +12,7 @@ couplingSystem/porosity/PIC.C
 couplingSystem/porosity/subDivision9.C
 couplingSystem/porosity/subDivision29.C
 couplingSystem/porosity/subDivision29Mod.C
+couplingSystem/porosity/diffusion.C
 couplingSystem/interaction/drag/drag.C
 couplingSystem/interaction/drag/ErgunWenYu.C
 couplingSystem/interaction/drag/DiFelice.C

phasicFlowCoupling/couplingSystem/couplingMesh.C

@@ -147,13 +147,13 @@ bool pFlow::coupling::couplingMesh::checkForDomainUpdate
         return true;
     }
 
-    if( abs(t-lastTimeUpdated_) < 0.98*fluidDt )
+    if( std::abs(t-lastTimeUpdated_) < static_cast<real>(0.98*fluidDt) )
     {
         lastTimeUpdated_ = t;
         return true;
     }
 
-    if( abs(t-(lastTimeUpdated_+domainUpdateInterval_)) < 0.98*fluidDt)
+    if( std::abs(t-(lastTimeUpdated_+domainUpdateInterval_)) < static_cast<real>(0.98*fluidDt))
     {
         lastTimeUpdated_ = t;
         return true;

phasicFlowCoupling/couplingSystem/porosity/PIC.C

@@ -48,24 +48,20 @@ bool pFlow::coupling::PIC::internalFieldUpdate()
 
 	auto& solidVol = solidVoldTmp.ref();
 
-	numInMesh_ = 0;
 	size_t numPar = centerMass_.size();
 
-#pragma omp parallel for reduction(+:numInMesh_)
+	#pragma omp parallel for
 	for(size_t i=0; i<numPar; i++)
 	{
-
-		auto cellId = cMesh_.findCellTree(centerMass_[i], parCellIndex_[i]);
+		const auto cellId = parCellIndex_[i];
 		if( cellId >= 0 )
 		{
 			#pragma omp atomic
 			solidVol[cellId] += 
 				static_cast<real>(3.14159265358979/6)*
 				pFlow::pow(particleDiameter_[i], static_cast<real>(3.0));
-				numInMesh_++;
+
 		}
-		parCellIndex_[i] = cellId;
-
 	}
 
 	this->ref() = Foam::max(

phasicFlowCoupling/couplingSystem/porosity/PIC.hpp

@@ -71,10 +71,6 @@ class PIC
 
 	bool internalFieldUpdate() override;
 
-	int32 numInMesh()const override
-	{
-		return numInMesh_;
-	}
 
 }; 
 

phasicFlowCoupling/couplingSystem/porosity/diffusion.C

@@ -0,0 +1,135 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+// from OpenFOAM
+#include "fvCFD.H"
+
+
+#include "diffusion.hpp"
+
+
+Foam::tmp<Foam::fvMatrix<Foam::scalar>> pFlow::coupling::diffusion::fvmDdt
+(
+    const Foam::volScalarField& sField
+)
+{
+	Foam::tmp<fvMatrix<Foam::scalar>> tfvm
+    (
+        new Foam::fvMatrix<Foam::scalar>
+        (
+            sField,
+            sField.dimensions()*Foam::dimVol/Foam::dimTime
+        )
+    );
+
+    Foam::fvMatrix<Foam::scalar>& fvm = tfvm.ref();
+
+    const Foam::scalar rDeltaT = 1.0/dt_.value();
+
+    fvm.diag() = rDeltaT*sField.mesh().Vsc();
+
+    if (sField.mesh().moving())
+    {
+        fvm.source() = rDeltaT*sField.oldTime().primitiveField()*sField.mesh().Vsc0();
+    }
+    else
+    {
+        fvm.source() = rDeltaT*sField.oldTime().primitiveField()*sField.mesh().Vsc();
+    }
+
+    return tfvm;
+}
+
+pFlow::coupling::diffusion::diffusion(
+	Foam::dictionary 		dict, 
+	couplingMesh& 			cMesh, 
+	MPI::centerMassField& 	centerMass, 
+	MPI::realProcCMField& 	parDiam)
+:
+	PIC(dict, cMesh, centerMass, parDiam),
+	nSteps_(Foam::max(1,dict.lookup<Foam::label>("nSteps"))),
+	intTime_(dict.lookup<Foam::scalar>("intTime")),
+	dt_("dt", Foam::dimTime, intTime_/nSteps_),
+	picSolDict_("picSolDict")
+{
+
+	picSolDict_.add("relTol", 0);
+	picSolDict_.add("tolerance", 1.0e-8);
+	picSolDict_.add("solver", "smoothSolver");
+	picSolDict_.add("smoother", "symGaussSeidel");
+}
+
+
+bool pFlow::coupling::diffusion::internalFieldUpdate()
+{
+
+	auto solidVoldTmp = Foam::volScalarField::Internal::New(
+		"solidVol",
+		this->mesh(),
+		Foam::dimensioned("solidVol", Foam::dimVol, Foam::scalar(0))
+	);
+
+	auto& solidVol = solidVoldTmp.ref();
+
+	size_t numPar = centerMass_.size();
+
+	#pragma omp parallel for
+	for(size_t i=0; i<numPar; i++)
+	{
+		const auto cellId = parCellIndex_[i];
+		if( cellId >= 0 )
+		{
+			#pragma omp atomic
+			solidVol[cellId] += 
+				static_cast<real>(3.14159265358979/6)*
+				pFlow::pow(particleDiameter_[i], static_cast<real>(3.0));
+
+		}
+	}
+
+	auto picAlphaTmp = Foam::volScalarField::New(
+		"picAlpha",
+		this->mesh(),
+		Foam::dimensioned("picAlpha", Foam::dimless, Foam::scalar(0)),
+		"zeroGradient"
+	);
+
+	volScalarField& picAlpha = picAlphaTmp.ref();
+
+
+	picAlpha.ref() = Foam::max( 1 - solidVol/this->mesh().V(), 0.0);
+	picAlpha.correctBoundaryConditions();
+
+
+	// start of Time loop
+	for(Foam::label i=0; i<nSteps_; i++)
+	{
+		picAlpha.storeOldTime();
+		fvScalarMatrix alphaEq
+		(
+			fvmDdt(picAlpha) - fvm::laplacian(picAlpha)
+		);
+		alphaEq.solve(picSolDict_);
+	}
+
+	this->ref() = picAlpha.internalField();
+
+	return true;
+}

phasicFlowCoupling/couplingSystem/porosity/diffusion.hpp

@@ -0,0 +1,91 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __diffusion_hpp__ 
+#define __diffusion_hpp__
+
+#include "virtualConstructor.hpp"
+
+// from phasicFlow-coupling
+#include "PIC.hpp"
+
+
+namespace pFlow::coupling
+{
+
+/**
+ * Particle In Cell (diffusion) model for porosity calculation
+ * 
+ * This model only considers the particle center and if the particle center 
+ * resides inside a cell, it is assumed that the whole volume of the particle
+ * is located in that cell.
+ * 
+ */
+class diffusion
+: 
+	public PIC
+{
+private:
+
+	Foam::label 			  nSteps_;
+
+	Foam::scalar 			  intTime_;
+
+	Foam::dimensionedScalar   dt_;
+
+	Foam::dictionary          picSolDict_;
+
+	Foam::tmp<Foam::fvMatrix<Foam::scalar>> fvmDdt
+    (
+        const Foam::volScalarField& sField
+    );
+
+public:
+
+	/// Type info
+	TypeInfo("diffusion");
+
+	/// Construc from dictionary 
+	diffusion(
+		Foam::dictionary 		dict, 
+		couplingMesh& 			cMesh, 
+		MPI::centerMassField& 	centerMass, 
+		MPI::realProcCMField& 	parDiam);
+
+	/// Destructor
+	virtual ~diffusion() = default;
+
+	/// Add this constructor to the list of virtual constructors
+	add_vCtor
+	(
+		porosity,
+		diffusion,
+		dictionary
+	);
+
+	bool internalFieldUpdate() override;
+
+
+}; 
+
+} // pFlow::coupling
+
+
+#endif

phasicFlowCoupling/couplingSystem/porosity/porosity.C

@@ -24,6 +24,23 @@ Licence:
 #include "streams.hpp"
 #include "Timer.hpp"
 
+
+void pFlow::coupling::porosity::mapCenters()
+{
+	numInMesh_ = 0;
+	size_t numPar = centerMass_.size();
+
+
+	#pragma omp parallel for reduction(+:numInMesh_)
+	for(size_t i = 0; i<numPar; i++)
+	{
+		auto cellId = cMesh_.findCellTree(centerMass_[i], parCellIndex_[i]);
+		parCellIndex_[i] = cellId;
+		if( cellId >= 0 ) numInMesh_++;	
+	}
+}
+
+
 pFlow::coupling::porosity::porosity(
 	Foam::dictionary		dict, 
 	couplingMesh& 			cMesh, 
@@ -61,6 +78,7 @@ void pFlow::coupling::porosity::calculatePorosity()
 {
 	Timer t;
 	t.start();
+	mapCenters();
 	this->internalFieldUpdate();
 	t.end();
 	output<<"mapping execution time " << t.lastTime()<<endl;

phasicFlowCoupling/couplingSystem/porosity/porosity.hpp

@@ -62,6 +62,10 @@ class porosity
 	/// cell indices of particles in this processor 
 	MPI::procCMField<Foam::label> 	parCellIndex_;
 
+	int32 										numInMesh_ = 0;
+
+	void mapCenters();
+
 public:
 
 	/// Type info
@@ -150,8 +154,10 @@ class porosity
 		bool internalFieldUpdate() = 0;
 
 		/// Return number of center mass points found in this mesh (processor)
-		virtual
-		int32 numInMesh()const = 0;
+		int32 numInMesh()const
+		{
+			return numInMesh_;
+		}
 
 		/// Report (output) number of center mass points found in all processors 
 		/// It is effective only in master processor