Merge pull request #61 from PhasicFlow/coupling

modification for coupling, data transfer

DEMSystems/CMakeLists.txt

@@ -2,6 +2,7 @@
 set(SourceFiles 
 DEMSystem/DEMSystem.cpp
 sphereDEMSystem/sphereDEMSystem.cpp
+sphereDEMSystem/sphereFluidParticles.cpp
 domainDistribute/domainDistribute.cpp
 )
 

DEMSystems/DEMSystem/DEMSystem.cpp

@@ -18,6 +18,9 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+// from phasicFlow
+#include "KokkosTypes.hpp"
+
 #include "DEMSystem.hpp"
 
 
@@ -27,10 +30,16 @@ pFlow::DEMSystem::DEMSystem(
 	int argc, 
 	char* argv[])
 :
-	ControlDict_(),
-	domains_(domains)
+	ControlDict_()
 {
   
+  REPORT(0)<<"Initializing host/device execution spaces . . . \n";
+	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
+	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
+
+ 	// initialize Kokkos
+	Kokkos::initialize( argc, argv ); 
+
 	REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
 	Control_ = makeUnique<systemControl>(
 		ControlDict_.startTime(),
@@ -43,7 +52,13 @@ pFlow::DEMSystem::DEMSystem(
 }
 
 pFlow::DEMSystem::~DEMSystem()
-{}
+{
+	
+	Control_.reset();
+
+	output<< "\nFinalizing host/device execution space ...."<<endl;
+	Kokkos::finalize();
+}
 
 
 pFlow::uniquePtr<pFlow::DEMSystem>

DEMSystems/DEMSystem/DEMSystem.hpp

@@ -43,16 +43,8 @@ protected:
 
 	uniquePtr<systemControl>		Control_ = nullptr;
 
-	std::vector<box> 				domains_;
-
 	uniquePtr<Timers> 				timers_;
 
-	// methods 
-	auto& Control()
-	{
-		return Control_();
-	}
-
 
 public:
 
@@ -91,6 +83,17 @@ public:
 		return Control_->g();
 	}
 
+	// methods 
+	auto& Control()
+	{
+		return Control_();
+	}
+
+	const auto& Control()const
+	{
+		return Control_();
+	}
+
 	auto inline constexpr usingDoulle()const
 	{
 		return pFlow::usingDouble__;
@@ -101,6 +104,9 @@ public:
 		return Control_->timers();
 	}
 
+	virtual
+	bool updateParticleDistribution(real extentFraction, const std::vector<box> domains) = 0;
+
 	virtual 
 	int32 numParInDomain(int32 di)const = 0;
 	
@@ -111,14 +117,32 @@ public:
 	span<const int32> parIndexInDomain(int32 di)const = 0;
 
 	virtual
-	bool changeDomainsSizeUpdateParticles(const std::vector<box>& domains) = 0;
+	span<real> parDiameter() = 0; 
 
 	virtual
-	bool updateParticles() = 0;
+	span<realx3> parVelocity() = 0;
+
+	virtual
+	span<realx3> parPosition() = 0;
+
+	virtual
+	span<realx3> parFluidForce() = 0;
+
+	virtual
+	span<realx3> parFluidTorque() = 0;
+
+	virtual 
+	bool sendFluidForceToDEM() = 0;
+
+	virtual
+	bool sendFluidTorqueToDEM() = 0;
 	
 	virtual 
 	real maxBounndingSphereSize()const = 0;
 
+	virtual 
+	bool beforeIteration() = 0;
+
 	virtual 
 	bool iterate(int32 n, real timeToWrite, word timeName) = 0;
 	

DEMSystems/domainDistribute/domainDistribute.cpp

@@ -49,6 +49,60 @@ maxBoundingBoxSize_(maxBoundingBox)
 bool pFlow::domainDistribute::locateParticles(
 		ViewType1D<realx3,HostSpace> points, includeMask mask)
 {
+	
+	range activeRange = mask.activeRange();
+
+	for(int32 di=0; di<numDomains_; di++)
+	{
+		particlesInDomains_[di].clear();
+	}
+
+	for(int32 i=activeRange.first; i<activeRange.second; i++)
+	{
+		if(mask(i))
+		{
+			for(int32 di=0; di<numDomains_; di++)
+			{
+				if(extDomains_[di].isInside(points[i]))
+				{
+					particlesInDomains_[di].push_back(i);
+				}
+			}
+		}
+	}
+
+	for(int32 di=0; di<numDomains_; di++)
+	{
+		numParInDomain_[di] = particlesInDomains_[di].size();
+	}
+
+	output<<" numParInDomain_ "<< numParInDomain_<<endl;
+
+	return true;
+}
+
+bool pFlow::domainDistribute::changeDomainsSize(
+		real extentFraction, 
+		real maxBoundingBoxSize, 
+		const std::vector<box>& domains)
+{
+	domainExtension_ = extentFraction;
+	maxBoundingBoxSize_ = maxBoundingBoxSize;
+
+	if(domains.size()!= numDomains_)
+	{
+		fatalErrorInFunction<<"number of new domians differs"<<endl;
+		return false;
+	}
+
+	clcDomains(domains);
+	return true;
+}
+
+
+/*bool pFlow::domainDistribute::locateParticles(
+		ViewType1D<realx3,HostSpace> points, includeMask mask)
+{
 	range active = mask.activeRange();
 	auto numInDomain = numParInDomain_.deviceVectorAll();
 	auto numDomains = numDomains_;
@@ -102,17 +156,4 @@ bool pFlow::domainDistribute::locateParticles(
 
 
 	return true;
-}
-
-bool pFlow::domainDistribute::changeDomainsSize(
-	const std::vector<box>& domains)
-{
-	if(domains.size()!= numDomains_)
-	{
-		fatalErrorInFunction<<"number of new domians differs"<<endl;
-		return false;
-	}
-
-	clcDomains(domains);
-	return true;
-}
+}*/

DEMSystems/domainDistribute/domainDistribute.hpp

@@ -83,7 +83,10 @@ public:
 			);
 	}
 
-	bool changeDomainsSize(const std::vector<box>& domains);
+	bool changeDomainsSize(
+		real extentFraction, 
+		real maxBoundingBoxSize, 
+		const std::vector<box>& domains);
 	
 
 	//template<typename includeMask>

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -30,15 +30,6 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 	DEMSystem(demSystemName, domains, argc, argv)
 { 
 	
-	REPORT(0)<<"Initializing host/device execution spaces . . . \n";
-	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
-	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
-
- 	// initialize Kokkos
-	Kokkos::initialize( argc, argv );  
-  
-	
-
 	REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
 	property_ = makeUnique<property>(
     	Control().caseSetup().path()+propertyFile__);
@@ -47,7 +38,7 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 	geometry_ = geometry::create(Control(), Property());
 
 	REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
-	particles_ = makeUnique<sphereParticles>(Control(), Property());
+	particles_ = makeUnique<sphereFluidParticles>(Control(), Property());
 
 
 	REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
@@ -65,18 +56,35 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 
 pFlow::sphereDEMSystem::~sphereDEMSystem()
 {
-	interaction_.reset();
-	insertion_.reset();
-	particles_.reset();
-	geometry_.reset();
-	property_.reset();
-	Control_.reset();
 	
-	output<< "\nFinalizing host/device execution space ...."<<endl;
-	Kokkos::finalize();
 }
 
+bool pFlow::sphereDEMSystem::updateParticleDistribution(
+	real extentFraction,
+	const std::vector<box> domains)
+{
 
+	if(!particleDistribution_->changeDomainsSize(
+		extentFraction,
+		maxBounndingSphereSize(),
+		domains))
+	{
+		fatalErrorInFunction<<
+		"cannot change the domain size"<<endl;
+		return false;
+	}
+
+	if(!particleDistribution_->locateParticles(
+			parPosition_, 
+			particles_->pStruct().activePointsMaskH()))
+	{
+		fatalErrorInFunction<<
+		"error in locating particles among sub-domains"<<endl;
+		return false;
+	}
+
+	return true;
+}
 
 pFlow::int32 
 	pFlow::sphereDEMSystem::numParInDomain(int32 di)const
@@ -104,24 +112,54 @@ pFlow::sphereDEMSystem::parIndexInDomain(int32 di)const
 	return particleDistribution_->particlesInDomain(di);
 }
 
-bool pFlow::sphereDEMSystem::changeDomainsSizeUpdateParticles(
-	const std::vector<box>& domains)
+pFlow::span<pFlow::real> pFlow::sphereDEMSystem::parDiameter() 
 {	
-	if( !particleDistribution_->changeDomainsSize(domains))
-		return false;
-
-	// should update list of particles here 
-	//************************************************************************************************
-	notImplementedFunction;
-	return false;	
+	return span<real>(parDiameter_.data(), parDiameter_.size());
 }
 
-bool pFlow::sphereDEMSystem::updateParticles()
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parVelocity()  
 {
-	notImplementedFunction;
-	return false;
+	return span<realx3>(parVelocity_.data(), parVelocity_.size());
 }
 
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parPosition()  
+{
+	return span<realx3>(parPosition_.data(), parPosition_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidForce() 
+{
+	auto& hVec =  particles_->fluidTorqueHostAll();
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::parFluidTorque() 
+{
+	auto& hVec =  particles_->fluidForceHostAll();
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+bool pFlow::sphereDEMSystem::sendFluidForceToDEM() 
+{
+	particles_->fluidForceHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::sphereDEMSystem::sendFluidTorqueToDEM()
+{
+	particles_->fluidTorqueHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::sphereDEMSystem::beforeIteration()
+{
+	parVelocity_ = particles_->dynPointStruct().velocityHostAll();
+	parPosition_ = particles_->dynPointStruct().pointPositionHostAll();
+	parDiameter_ = particles_->diameter().hostVectorAll();
+	return true;
+}
+
+
 bool pFlow::sphereDEMSystem::iterate(
 	int32 n, 
 	real timeToWrite, 

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -25,7 +25,7 @@ Licence:
 #include "property.hpp"
 #include "uniquePtr.hpp"
 #include "geometry.hpp"
-#include "sphereParticles.hpp"
+#include "sphereFluidParticles.hpp"
 #include "interaction.hpp"
 #include "Insertions.hpp"
 #include "domainDistribute.hpp"
@@ -46,7 +46,7 @@ protected:
 
 	uniquePtr<geometry> 		geometry_ = nullptr;
 
-	uniquePtr<sphereParticles> 	particles_ = nullptr;
+	uniquePtr<sphereFluidParticles> 	particles_ = nullptr;
 
 	uniquePtr<sphereInsertion> 	insertion_ = nullptr;
 
@@ -54,6 +54,13 @@ protected:
 
 	uniquePtr<domainDistribute> particleDistribution_=nullptr;
 
+	// to be used for CPU communications 
+	ViewType1D<realx3, HostSpace>  parVelocity_;
+
+	ViewType1D<realx3, HostSpace>  parPosition_;
+
+	ViewType1D<real, HostSpace> parDiameter_;
+
 //  protected member functions
 	auto& Property()
 	{
@@ -97,6 +104,7 @@ public:
 		word);
 
 
+	bool updateParticleDistribution(real extentFraction, const std::vector<box> domains) override;
 	 
 	int32 numParInDomain(int32 di)const override;
 	
@@ -104,9 +112,21 @@ public:
 
 	span<const int32> parIndexInDomain(int32 di)const override;
 
-	bool changeDomainsSizeUpdateParticles(const std::vector<box>& domains) override;
+	span<real> parDiameter() override;
 
-	bool updateParticles() override;	
+	span<realx3> parVelocity() override;
+
+	span<realx3> parPosition() override;
+
+	span<realx3> parFluidForce() override;
+
+	span<realx3> parFluidTorque() override;
+
+	bool sendFluidForceToDEM() override;
+
+	bool sendFluidTorqueToDEM() override;
+
+	bool beforeIteration() override;
 
 	bool iterate(int32 n, real timeToWrite, word timeName) override;
 	

DEMSystems/sphereDEMSystem/sphereFluidParticles.cpp

@@ -0,0 +1,100 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+
+#include "sphereFluidParticles.hpp"
+#include "sphereFluidParticlesKernels.hpp"
+
+
+pFlow::sphereFluidParticles::sphereFluidParticles(
+	systemControl &control,
+	const property& prop
+)
+:
+	sphereParticles(control, prop),
+	fluidForce_(
+		this->time().emplaceObject<realx3PointField_HD>(
+			objectFile(
+				"fluidForce",
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_ALWAYS
+				),
+			pStruct(),
+			zero3
+			)
+		),
+	fluidTorque_(
+		this->time().emplaceObject<realx3PointField_HD>(
+			objectFile(
+				"fluidTorque",
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_ALWAYS
+				),
+			pStruct(),
+			zero3
+			)
+		)
+{} 
+
+
+
+bool pFlow::sphereFluidParticles::iterate() 
+{
+
+	accelerationTimer_.start();
+		pFlow::sphereFluidParticlesKernels::acceleration(
+			control().g(),
+			mass().deviceVectorAll(),
+			contactForce().deviceVectorAll(),
+			fluidForce().deviceVectorAll(),
+			I().deviceVectorAll(),
+			contactTorque().deviceVectorAll(),
+			fluidTorque().deviceVectorAll(),
+			pStruct().activePointsMaskD(),
+			accelertion().deviceVectorAll(),
+			rAcceleration().deviceVectorAll()
+			);
+	accelerationTimer_.end();
+	
+	intCorrectTimer_.start();
+	
+		dynPointStruct_.correct(this->dt(), accelertion_);
+	
+		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
+	
+	intCorrectTimer_.end();
+	
+	return true;
+}
+
+void pFlow::sphereFluidParticles::fluidForceHostUpdatedSync()
+{
+	fluidForce_.modifyOnHost();
+	fluidForce_.syncViews();
+	return;
+}
+
+void pFlow::sphereFluidParticles::fluidTorqueHostUpdatedSync()
+{
+	fluidTorque_.modifyOnHost();
+	fluidTorque_.syncViews();
+	return;
+}

DEMSystems/sphereDEMSystem/sphereFluidParticles.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.
+
+-----------------------------------------------------------------------------*/
+
+/*! 
+   @class pFlow::sphereParticles
+    
+   @brief Class for managing spherical particles 
+   
+   This is a top-level class that contains the essential components for
+   defining spherical prticles in a DEM simulation.
+ */
+
+#ifndef __sphereFluidParticles_hpp__
+#define __sphereFluidParticles_hpp__ 
+
+#include "sphereParticles.hpp"
+
+
+namespace pFlow
+{
+
+class sphereFluidParticles
+:
+	public sphereParticles 
+{
+protected:
+	
+	/// pointField of rotational acceleration of particles on device  
+	realx3PointField_HD& 	fluidForce_;
+
+	realx3PointField_HD& 	fluidTorque_;
+
+
+public:
+
+	/// construct from systemControl and property 
+	sphereFluidParticles(systemControl &control, const property& prop);
+	
+
+	/// iterate particles 
+	bool iterate() override;	
+
+	
+	auto& fluidForce()
+	{
+		return fluidForce_;
+	}
+
+	auto& fluidTorque()
+	{
+		return fluidTorque_;
+	}
+
+
+	auto& fluidForceHostAll()
+	{
+		return fluidForce_.hostVectorAll();
+	}
+
+	auto& fluidTorqueHostAll()
+	{
+		return fluidTorque_.hostVectorAll();
+	}
+
+
+	void fluidForceHostUpdatedSync();
+
+	void fluidTorqueHostUpdatedSync();
+	
+}; //sphereFluidParticles
+
+} // pFlow
+
+#endif //__sphereFluidParticles_hpp__

DEMSystems/sphereDEMSystem/sphereFluidParticlesKernels.hpp

@@ -0,0 +1,79 @@
+/*------------------------------- 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 __sphereFluidParticlesKernels_hpp__
+#define __sphereFluidParticlesKernels_hpp__
+
+namespace pFlow::sphereFluidParticlesKernels
+{
+
+using rpAcceleration = Kokkos::RangePolicy<
+			DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<int32>
+			>;
+
+template<typename IncludeFunctionType>
+void acceleration( 
+	realx3		g,
+	deviceViewType1D<real>  	mass,
+	deviceViewType1D<realx3>  	force,
+	deviceViewType1D<realx3>  	fluidForce,
+	deviceViewType1D<real>  	I,
+	deviceViewType1D<realx3>  	torque,
+	deviceViewType1D<realx3>  	fluidTorque,
+	IncludeFunctionType 		incld,
+	deviceViewType1D<realx3> 	lAcc,
+	deviceViewType1D<realx3> 	rAcc
+	)
+{
+
+	auto activeRange = incld.activeRange();
+	if(incld.allActive())
+	{
+		Kokkos::parallel_for(
+		"pFlow::sphereParticlesKernels::acceleration",
+		rpAcceleration(activeRange.first, activeRange.second),
+		LAMBDA_HD(int32 i){
+				lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
+				rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
+		});
+	}
+	else
+	{
+		Kokkos::parallel_for(
+		"pFlow::sphereParticlesKernels::acceleration",
+		rpAcceleration(activeRange.first, activeRange.second),
+		LAMBDA_HD(int32 i){
+			if(incld(i))
+			{
+				lAcc[i] = (force[i]+fluidForce[i])/mass[i] + g;
+				rAcc[i] = (torque[i]+fluidTorque[i])/I[i];
+			}
+		});
+
+	}
+	
+	Kokkos::fence();
+}
+
+}
+
+#endif 

src/Particles/dynamicPointStructure/dynamicPointStructure.hpp

@@ -100,6 +100,16 @@ public:
 		return velocity_;
 	}
 
+	inline auto velocityHostAll()
+	{
+		return velocity_.hostVectorAll();
+	}
+
+	inline auto pointPositionHostAll()
+	{
+		return pStruct_.pointPositionHostAll();
+	}
+
 	auto markDeleteOutOfBox(const box& domain)
 	{
 		return pStruct_.markDeleteOutOfBox(domain);

src/phasicFlow/repository/Time/timeControl.hpp

@@ -108,6 +108,10 @@ public:
 		return tmp;
 	}
 
+	real startTime()const
+	{
+		return startTime_;
+	}
 
 	real currentTime() const 
 	{

src/phasicFlow/structuredData/pointStructure/pointStructure.hpp

@@ -255,6 +255,12 @@ public:
 		FUNCTION_H
 		const int8Field_HD& pointFlag()const;
 
+		INLINE_FUNCTION_H
+		auto pointPositionHostAll()
+		{
+			return pointPosition_.hostVectorAll();
+		}
+
 		// - size of data structure
 		FUNCTION_H
 		label size()const;