grainDEMSystem folder is added

DEMSystems/CMakeLists.txt

@@ -4,8 +4,8 @@ domainDistribute/domainDistribute.cpp
 DEMSystem/DEMSystem.cpp
 sphereDEMSystem/sphereFluidParticles.cpp
 sphereDEMSystem/sphereDEMSystem.cpp
-#
-
+grainDEMSystem/grainFluidParticles.cpp
+grainDEMSystem/grainDEMSystem.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)

DEMSystems/DEMSystem/DEMSystem.hpp

@@ -121,6 +121,9 @@ public:
 	virtual
 	span<real> diameter() = 0; 
 
+	virtual 
+	span<real> courseGrainFactor() = 0;
+
 	virtual
 	span<realx3> acceleration()=0;
 

DEMSystems/grainDEMSystem/grainDEMSystem.cpp

@@ -0,0 +1,268 @@
+/*------------------------------- 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 "grainDEMSystem.hpp"
+#include "vocabs.hpp"
+
+bool pFlow::grainDEMSystem::loop()
+{
+
+	do 
+	{
+		//
+		if(! insertion_().insertParticles( 
+			Control().time().currentIter(),
+			Control().time().currentTime(),
+			Control().time().dt()	) )
+		{
+			fatalError<<
+			"particle insertion failed in grainDEMSystem.\n";
+			return false;
+		}	
+
+		geometry_->beforeIteration();
+		
+		interaction_->beforeIteration();
+
+		particles_->beforeIteration();
+		
+		
+		interaction_->iterate();
+
+		particles_->iterate();
+
+		geometry_->iterate();
+
+		particles_->afterIteration();
+		
+		geometry_->afterIteration();
+
+	}while(Control()++);
+
+	return true;
+}
+
+pFlow::grainDEMSystem::grainDEMSystem(
+		word  demSystemName,
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[],
+		bool requireRVel)
+:
+	DEMSystem(demSystemName, domains, argc, argv),
+	requireRVel_(requireRVel)
+{ 
+	
+	REPORT(0)<<"\nReading proprties . . . "<<END_REPORT;
+	property_ = makeUnique<property>(
+		propertyFile__,
+    	Control().caseSetup().path());
+
+	REPORT(0)<< "\nCreating surface geometry for grainDEMSystem . . . "<<END_REPORT;
+	geometry_ = geometry::create(Control(), Property());
+
+
+	REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
+	particles_ = makeUnique<grainFluidParticles>(Control(), Property());
+
+
+	insertion_ = makeUnique<grainInsertion>( 
+		particles_(), 
+		particles_().grains());
+	
+	REPORT(0)<<"\nCreating interaction model for grain-grain contact and grain-wall contact . . ."<<END_REPORT;
+	interaction_ = interaction::create(
+		Control(),
+		Particles(),
+		Geometry());
+
+	auto maxD = maxBounndingSphereSize();
+
+	particleDistribution_ = makeUnique<domainDistribute>(domains, maxD);
+
+}
+
+
+pFlow::grainDEMSystem::~grainDEMSystem()
+{
+	
+}
+
+bool pFlow::grainDEMSystem::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(
+			positionHost_, 
+			particles_->pStruct().activePointsMaskHost()))
+	{
+		fatalErrorInFunction<<
+		"error in locating particles among sub-domains"<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+pFlow::int32 
+	pFlow::grainDEMSystem::numParInDomain(int32 di)const
+{
+	return particleDistribution_().numParInDomain(di);
+}
+
+std::vector<pFlow::int32> 
+	pFlow::grainDEMSystem::numParInDomains()const
+{
+	const auto& distribute = particleDistribution_();
+	int32 numDomains = distribute.numDomains();
+	std::vector<int32> nums(numDomains);
+	for(int32 i=0; i<numDomains; i++)
+	{
+		nums[i] = distribute.numParInDomain(i);
+	}
+
+	return nums;
+}
+
+pFlow::span<const pFlow::int32> 
+pFlow::grainDEMSystem::parIndexInDomain(int32 di)const
+{
+	return particleDistribution_->particlesInDomain(di);
+}
+
+pFlow::span<pFlow::real> pFlow::grainDEMSystem::diameter() 
+{	
+	return span<real>(diameterHost_.data(), diameterHost_.size());
+}
+
+pFlow::span<pFlow::real> pFlow::grainDEMSystem::courseGrainFactor()
+{
+    return span<real>(particles_->courseGrainFactorHost().data(), particles_->courseGrainFactorHost().size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::acceleration()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::velocity()  
+{
+	return span<realx3>(velocityHost_.data(), velocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::position()  
+{
+	return span<realx3>(positionHost_.data(), positionHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rAcceleration()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rVelocity()
+{
+    return span<realx3>(rVelocityHost_.data(), rVelocityHost_.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::rPosition()
+{
+    return span<realx3>(nullptr, 0);
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidForce() 
+{
+	auto& hVec =  particles_->fluidForceHost();
+	
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+pFlow::span<pFlow::realx3> pFlow::grainDEMSystem::parFluidTorque() 
+{
+	auto& hVec =  particles_->fluidTorqueHost();
+	return span<realx3>(hVec.data(), hVec.size());
+}
+
+bool pFlow::grainDEMSystem::sendFluidForceToDEM() 
+{
+	particles_->fluidForceHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::grainDEMSystem::sendFluidTorqueToDEM()
+{
+	particles_->fluidTorqueHostUpdatedSync();
+	return true;
+}
+
+bool pFlow::grainDEMSystem::beforeIteration()
+{
+	velocityHost_ = std::as_const(particles_()).velocity().hostView();
+	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
+	diameterHost_ = particles_->diameter().hostView();
+
+	if(requireRVel_)
+		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
+	
+
+	return true;
+}
+
+
+bool pFlow::grainDEMSystem::iterate(
+	real upToTime, 
+	real timeToWrite, 
+	word timeName) 
+{
+
+	Control().time().setStopAt(upToTime);
+	Control().time().setOutputToFile(timeToWrite, timeName);
+	
+	return loop();
+
+	return true;
+}
+
+bool pFlow::grainDEMSystem::iterate(real upToTime)
+{
+	Control().time().setStopAt(upToTime);
+	return loop();
+	return true;
+}
+
+pFlow::real 
+	pFlow::grainDEMSystem::maxBounndingSphereSize()const
+{
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	return maxD;
+}

DEMSystems/grainDEMSystem/grainDEMSystem.hpp

@@ -0,0 +1,163 @@
+/*------------------------------- 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 __grainDEMSystem_hpp__
+#define __grainDEMSystem_hpp__
+
+#include "DEMSystem.hpp"
+#include "property.hpp"
+#include "uniquePtr.hpp"
+#include "geometry.hpp"
+#include "grainFluidParticles.hpp"
+#include "interaction.hpp"
+#include "Insertions.hpp"
+#include "domainDistribute.hpp"
+
+
+namespace pFlow
+{
+
+class grainDEMSystem
+:
+	public DEMSystem
+{
+protected:
+	
+	// protected members 
+
+	uniquePtr<property> 			property_ = nullptr;
+
+	uniquePtr<geometry> 			geometry_ = nullptr;
+
+	uniquePtr<grainFluidParticles>  particles_ = nullptr;
+
+	uniquePtr<grainInsertion> 		insertion_ = nullptr;
+
+	uniquePtr<interaction> 			interaction_ = nullptr;
+
+	uniquePtr<domainDistribute> 	particleDistribution_=nullptr;
+
+	// to be used for CPU communications 
+	ViewType1D<realx3, HostSpace> 	velocityHost_;
+
+	ViewType1D<realx3, HostSpace> 	positionHost_;
+
+	ViewType1D<real, HostSpace> 	diameterHost_;
+
+	bool 							requireRVel_ = false;
+
+	ViewType1D<realx3, HostSpace> 	rVelocityHost_;
+
+
+//  protected member functions
+	auto& Property()
+	{
+		return property_();
+	}
+
+	auto& Geometry()
+	{
+		return geometry_();
+	}
+
+	auto& Particles()
+	{
+		return particles_();
+	}
+
+	auto& Interaction()
+	{
+		return interaction_();
+	} 
+
+	bool loop();
+
+public:
+
+	TypeInfo("grainDEMSystem");
+
+	grainDEMSystem(
+		word  demSystemName,
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[],
+		bool requireRVel=false);
+
+	virtual ~grainDEMSystem();
+
+	grainDEMSystem(const grainDEMSystem&)=delete;
+
+	grainDEMSystem& operator = (const grainDEMSystem&)=delete;
+
+	add_vCtor(
+		DEMSystem,
+		grainDEMSystem,
+		word);
+
+
+	bool updateParticleDistribution(real extentFraction, const std::vector<box> domains) override;
+	 
+	int32 numParInDomain(int32 di)const override;
+	
+	std::vector<int32> numParInDomains()const override;
+
+	span<const int32> parIndexInDomain(int32 di)const override;
+
+	span<real> diameter() override;
+
+	span<real> courseGrainFactor() override;
+
+	span<realx3> acceleration() override;
+
+	span<realx3> velocity() override;
+
+	span<realx3> position() override;
+
+	span<realx3> rAcceleration() override;
+
+	span<realx3> rVelocity() override;
+
+	span<realx3> rPosition() override;
+
+	span<realx3> parFluidForce() override;
+
+	span<realx3> parFluidTorque() override;
+
+	bool sendFluidForceToDEM() override;
+
+	bool sendFluidTorqueToDEM() override;
+
+	bool beforeIteration() override;
+
+	bool iterate(
+		real upToTime, 
+		real timeToWrite, 
+		word timeName) override;
+
+	bool iterate(real upToTime) override;
+	
+	real maxBounndingSphereSize()const override;
+
+};
+
+
+} // pFlow
+
+#endif

DEMSystems/grainDEMSystem/grainFluidParticles.cpp

@@ -0,0 +1,103 @@
+/*------------------------------- 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 "grainFluidParticles.hpp"
+#include "grainFluidParticlesKernels.hpp"
+
+
+pFlow::grainFluidParticles::grainFluidParticles(
+	systemControl &control,
+	const property& prop
+)
+:
+	grainParticles(control, prop),
+	fluidForce_(
+    	objectFile(
+        	"fluidForce",
+    		"",
+    		objectFile::READ_IF_PRESENT,
+        	objectFile::WRITE_ALWAYS
+      	),
+      dynPointStruct(),
+      zero3
+    ),
+	fluidTorque_(
+		objectFile(
+			"fluidTorque",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
+      	),
+    	dynPointStruct(),
+    	zero3
+	)
+{} 
+
+bool pFlow::grainFluidParticles::beforeIteration()
+{
+	grainParticles::beforeIteration();
+	
+	// copy the data (if required) from device to host
+	courseGrainFactorHost_ = coarseGrainFactor().hostView();
+
+	return true;
+}
+
+
+bool pFlow::grainFluidParticles::iterate() 
+{
+
+	accelerationTimer().start();
+		pFlow::grainFluidParticlesKernels::acceleration(
+			control().g(),
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			fluidForce_.deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			fluidTorque_.deviceViewAll(),
+			pStruct().activePointsMaskDevice(),
+			accelertion().deviceViewAll(),
+			rAcceleration().deviceViewAll()
+			);
+	accelerationTimer().end();
+	
+	intCorrectTimer().start();
+	
+		dynPointStruct().correct(this->dt(), accelertion());
+	
+		rVelIntegration().correct(this->dt(), rVelocity(), rAcceleration());
+	
+	intCorrectTimer().end();
+	
+	return true;
+}
+
+void pFlow::grainFluidParticles::fluidForceHostUpdatedSync()
+{
+	copy(fluidForce_.deviceView(), fluidForceHost_);
+	return;
+}
+
+void pFlow::grainFluidParticles::fluidTorqueHostUpdatedSync()
+{
+	copy(fluidTorque_.deviceView(), fluidTorqueHost_);
+	return;
+}

DEMSystems/grainDEMSystem/grainFluidParticles.hpp

@@ -0,0 +1,102 @@
+/*------------------------------- 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::grainFluidParticles
+    
+   @brief Class for managing grain particles with fluid interactions 
+   
+   This is a top-level class that contains the essential components for
+   defining grain prticles with fluid interactions in a DEM simulation.
+ */
+
+#ifndef __grainFluidParticles_hpp__
+#define __grainFluidParticles_hpp__ 
+
+#include "grainParticles.hpp"
+
+
+namespace pFlow
+{
+
+class grainFluidParticles
+:
+	public grainParticles 
+{
+protected:
+	
+	/// pointField of rotational acceleration of particles on device  
+	realx3PointField_D 		fluidForce_;
+
+	hostViewType1D<realx3> 	fluidForceHost_;
+
+	realx3PointField_D 		fluidTorque_;
+
+	hostViewType1D<realx3>	fluidTorqueHost_;
+
+	hostViewType1D<real> 	courseGrainFactorHost_;
+
+public:
+
+	/// construct from systemControl and property 
+	grainFluidParticles(systemControl &control, const property& prop);
+
+	~grainFluidParticles() override = default;
+	
+	/// before iteration step 
+	bool beforeIteration() override;
+	
+	/// iterate particles 
+	bool iterate() override;	
+
+	auto& fluidForce()
+	{
+		return fluidForce_;
+	}
+
+	auto& fluidTorque()
+	{
+		return fluidTorque_;
+	}
+
+	auto& fluidForceHost()
+	{
+		return fluidForceHost_; 
+	}
+
+	auto& fluidTorqueHost()
+	{
+		return fluidTorqueHost_;
+	}
+
+	auto& courseGrainFactorHost()
+	{
+		return courseGrainFactorHost_;
+	}
+
+	void fluidForceHostUpdatedSync();
+
+	void fluidTorqueHostUpdatedSync();
+	
+}; //grainFluidParticles
+
+} // pFlow
+
+#endif //__sphereFluidParticles_hpp__

DEMSystems/grainDEMSystem/grainFluidParticlesKernels.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 __grainFluidParticlesKernels_hpp__
+#define __grainFluidParticlesKernels_hpp__
+
+namespace pFlow::grainFluidParticlesKernels
+{
+
+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.isAllActive())
+	{
+		Kokkos::parallel_for(
+		"pFlow::grainParticlesKernels::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::grainParticlesKernels::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 

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -164,6 +164,11 @@ pFlow::span<pFlow::real> pFlow::sphereDEMSystem::diameter()
 	return span<real>(diameterHost_.data(), diameterHost_.size());
 }
 
+pFlow::span<pFlow::real> pFlow::sphereDEMSystem::courseGrainFactor()
+{
+    return span<real>(nullptr, 0);
+}
+
 pFlow::span<pFlow::realx3> pFlow::sphereDEMSystem::acceleration()
 {
     return span<realx3>(nullptr, 0);

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -122,6 +122,8 @@ public:
 
 	span<real> diameter() override;
 
+	span<real> courseGrainFactor() override;
+
 	span<realx3> acceleration() override;
 
 	span<realx3> velocity() override;

src/Particles/GrainParticles/grainParticles/grainParticles.hpp

@@ -100,15 +100,23 @@ private:
 	  realVector&       Is,
 	  uint32Vector&     shIndex
 	);
-	/*bool initializeParticles();
+	
+protected:
 
-	bool insertSphereParticles(
-	    const wordVector& names,
-	    const int32IndexContainer& indices,
-	    bool setId = true);
+	Timer& accelerationTimer()
+	{
+		return accelerationTimer_;
+	}
 
-	virtual uniquePtr<List<eventObserver*>> getFieldObjectList()const override;
-	*/
+	Timer& intCorrectTimer()
+	{
+		return intCorrectTimer_;
+	}
+
+	integration& rVelIntegration()
+	{
+		return rVelIntegration_();
+	}
 
 public:
 
@@ -185,7 +193,12 @@ public:
 		return grainDiameter_;
 	}
 
-		const realPointField_D& coarseGrainFactor() const 
+	const realPointField_D& grainDiameter()const
+	{
+		return grainDiameter_;
+	}
+
+	const realPointField_D& coarseGrainFactor() const 
 	{
 		return coarseGrainFactor_;
 	}