sphereParticles tested

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -21,9 +21,10 @@ Licence:
 #include "sphereParticles.hpp"
 #include "systemControl.hpp"
 #include "vocabs.hpp"
-//#include "setFieldList.hpp"
-//#include "sphereParticlesKernels.hpp"
+#include "sphereParticlesKernels.hpp"
 
+
+//#include "setFieldList.hpp"
 /*pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
 pFlow::sphereParticles::getFieldObjectList()const
 {
@@ -205,16 +206,36 @@ bool pFlow::sphereParticles::initInertia()
 			exeSpace,
 			Kokkos::IndexType<uint32>>;
 	
-	auto aPointsMask = dynPointStruct().activePointsMaskDevice();
-	auto aRange = aPointsMask.activeRange();
+	auto [minIndex, maxIndex] = minMax(shapeIndex().internal());		
 
-	auto field_I = I_.fieldDevice();
-	auto field_shapeIndex = shapeIndex().fieldDevice();
-	
-	const auto& shp = getShapes();
-	realVector_D I ("I", shp.Inertia());
-	auto d_I = I.deviceVector();
+	if( !spheres_.indexValid(maxIndex) )
+	{
+		fatalErrorInFunction<< 
+		"the maximum value of shapeIndex is "<< maxIndex << 
+		" which is not valid."<<endl;
+		return false;
+	}
+
+	auto aPointsMask 	= dynPointStruct().activePointsMaskDevice();
+	auto aRange 		= aPointsMask.activeRange();
 	
+	auto field_shapeIndex 	= shapeIndex().fieldDevice();
+	auto field_diameter 	= diameter_.fieldDevice();
+	auto field_mass 		= mass_.fieldDevice();
+	auto field_propId 		= propertyId_.fieldDevice();
+	auto field_I 			= I_.fieldDevice();
+
+	// get info from spheres shape 
+	realVector_D d("diameter", spheres_.boundingDiameter());
+	realVector_D mass("mass",spheres_.mass());
+	uint32Vector_D propId("propId", spheres_.shapePropertyIds());
+	realVector_D  I("I", spheres_.Inertia());
+
+	auto d_d 		= d.deviceVector();
+	auto d_mass 	= mass.deviceVector();
+	auto d_propId 	= propId.deviceVector();
+	auto d_I 		= I.deviceVector();
+		
 	Kokkos::parallel_for(
 		"particles::initInertia",
 		policy(aRange.start(), aRange.end()),
@@ -224,9 +245,13 @@ bool pFlow::sphereParticles::initInertia()
 			{
 				uint32 index = field_shapeIndex[i];
 				field_I[i] = d_I[index];
+				field_diameter[i] = d_d[index];
+				field_mass[i] = d_mass[index];
+				field_propId[i] = d_propId[index];
 			}
 		});
 	Kokkos::fence();
+	
 	return true;
 }
 
@@ -243,6 +268,38 @@ pFlow::sphereParticles::sphereParticles(
 		&control.caseSetup(),
 		prop
 	),
+	propertyId_
+	(
+		objectFile
+		(
+			"propertyId",
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_NEVER
+		),
+		dynPointStruct(),
+		0u
+	),
+	diameter_
+	(
+		objectFile(
+			"diameter",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER),
+		dynPointStruct(),
+		0.00000000001
+	),
+	mass_
+	(
+		objectFile(
+			"mass",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER),
+		dynPointStruct(),
+		0.0000000001
+	),
 	I_
 	(
 		objectFile
@@ -285,7 +342,7 @@ pFlow::sphereParticles::sphereParticles(
 	intCorrectTimer_(
 		"Integration-correct", &this->timers() )
 {
-
+	
 	auto intMethod = control.settingsDict().getVal<word>("integrationMethod");
 	REPORT(1)<<"Creating integration method "<<Green_Text(intMethod)
 		  << " for rotational velocity."<<END_REPORT;
@@ -306,31 +363,10 @@ pFlow::sphereParticles::sphereParticles(
 	}
 
 	WARNING<<"setFields for rVelIntegration_"<<END_WARNING;
-	/*if(rVelIntegration_->needSetInitialVals())
-	{
-		
-		auto [minInd, maxInd] = pStruct().activeRange();
-		int32IndexContainer indexHD(minInd, maxInd);
-	
-		auto n = indexHD.size();
-		auto index = indexHD.indicesHost();
-
-		realx3Vector rvel(n,RESERVE());
-		const auto hrVel = rVelocity_.hostVector();
-
-		for(auto i=0; i<n; i++)
-		{
-			rvel.push_back( hrVel[index(i)]);
-		}
-		
-		REPORT(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endREPORT;
-		rVelIntegration_->setInitialVals(indexHD, rvel);
-
-	}*/
-
 	
 	if(!initInertia())
 	{
+		fatalErrorInFunction;
 		fatalExit;
 	}
 	
@@ -437,18 +473,16 @@ bool pFlow::sphereParticles::iterate()
 
 	particles::iterate();
 	accelerationTimer_.start();
-	WARNING<<"pFlow::sphereParticlesKernels::acceleration"<<END_WARNING;
-	//INFO<<"before accelerationTimer_ "<<endINFO;
-		/*pFlow::sphereParticlesKernels::acceleration(
+		pFlow::sphereParticlesKernels::acceleration(
 			control().g(),
-			mass().deviceVectorAll(),
-			contactForce().deviceVectorAll(),
-			I().deviceVectorAll(),
-			contactTorque().deviceVectorAll(),
-			pStruct().activePointsMaskD(),
-			accelertion().deviceVectorAll(),
-			rAcceleration().deviceVectorAll()
-			);*/
+			mass().fieldDevice(),
+			contactForce().fieldDevice(),
+			I().fieldDevice(),
+			contactTorque().fieldDevice(),
+			dynPointStruct().activePointsMaskDevice(),
+			accelertion().fieldDevice(),
+			rAcceleration().fieldDevice()
+			);
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
@@ -485,3 +519,13 @@ const pFlow::shape &pFlow::sphereParticles::getShapes() const
 {
     return spheres_;
 }
+
+void pFlow::sphereParticles::boundingSphereMinMax
+(
+	real & minDiam, 
+	real& maxDiam
+)const
+{
+	minDiam = spheres_.minBoundingSphere();
+	maxDiam = spheres_.maxBoundingSphere();
+}

src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp

@@ -51,6 +51,15 @@ private:
 		/// reference to shapes 
 	ShapeType 			spheres_;
 
+	/// property id on device
+	uint32PointField_D 		propertyId_;
+
+	/// diameter / boundig sphere size of particles on device
+	realPointField_D		diameter_;
+
+	/// mass of particles field 
+	realPointField_D 		mass_;
+
 	/// pointField of inertial of particles 
 	realPointField_D     	I_;
 
@@ -150,6 +159,24 @@ public:
 		notImplementedFunction;
 		return false;
 	}
+
+	
+	const uint32PointField_D& propertyId()const override
+	{
+		return propertyId_;
+	}
+	
+	
+	const realPointField_D& diameter()const override
+	{
+		return diameter_;
+	}
+
+
+	const realPointField_D& mass()const override
+	{
+		return mass_;
+	}
 	
 	/// before iteration step 
 	bool beforeIteration() override;
@@ -178,6 +205,9 @@ public:
 	word shapeTypeName()const override;
 	
 	const shape& getShapes()const override;
+
+	
+	void boundingSphereMinMax(real & minDiam, real& maxDiam)const override;
 	
 }; //sphereParticles
 

src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.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 "sphereParticlesKernels.hpp"
+
+using policy = Kokkos::RangePolicy<
+			pFlow::DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<pFlow::uint32>>;
+
+void pFlow::sphereParticlesKernels::addMassDiamInertiaProp
+(
+    deviceViewType1D<uint32>    shapeIndex,
+    deviceViewType1D<real>  	mass,
+	deviceViewType1D<real>  	diameter,
+	deviceViewType1D<real>  	I,
+    deviceViewType1D<uint32>  	propertyId,
+    pFlagTypeDevice 		    incld,
+    deviceViewType1D<real>  	src_mass,
+	deviceViewType1D<real>  	src_diameter,
+	deviceViewType1D<real>  	src_I,
+    deviceViewType1D<uint32>  	src_propertyId
+)
+{
+    auto aRange = incld.activeRange();
+
+    Kokkos::parallel_for(
+		"particles::initInertia",
+		policy(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(incld(i))
+			{
+				uint32 index = shapeIndex[i];
+				I[i] = src_I[index];
+				diameter[i] = src_diameter[index];
+				mass[i] = src_mass[index];
+				propertyId[i] = src_propertyId[index];
+			}
+		});
+
+}
+
+void pFlow::sphereParticlesKernels::acceleration
+( 
+	realx3		g,
+	deviceViewType1D<real>  	mass,
+	deviceViewType1D<realx3>  	force,
+	deviceViewType1D<real>  	I,
+	deviceViewType1D<realx3>  	torque,
+	pFlagTypeDevice 		    incld,
+	deviceViewType1D<realx3> 	lAcc,
+	deviceViewType1D<realx3> 	rAcc
+)
+{
+
+	auto activeRange = incld.activeRange();
+	if(incld.isAllActive())
+	{
+		Kokkos::parallel_for(
+		"pFlow::sphereParticlesKernels::acceleration",
+		policy(activeRange.start(), activeRange.end()),
+		LAMBDA_HD(uint32 i){
+				lAcc[i] = force[i]/mass[i] + g;
+				rAcc[i] = torque[i]/I[i];
+		});
+	}
+	else
+	{
+		Kokkos::parallel_for(
+		"pFlow::sphereParticlesKernels::acceleration",
+		policy(activeRange.start(), activeRange.end()),
+		LAMBDA_HD(uint32 i){
+			if(incld(i))
+			{
+				lAcc[i] = force[i]/mass[i] + g;
+				rAcc[i] = torque[i]/I[i];
+			}
+		});
+
+	}
+	
+	Kokkos::fence();
+}

src/Particles/SphereParticles/sphereParticles/sphereParticlesKernels.hpp

@@ -21,56 +21,36 @@ Licence:
 #ifndef __sphereParticlesKernels_hpp__
 #define __sphereParticlesKernels_hpp__
 
+#include "types.hpp"
+#include "pointFlag.hpp"
+
 namespace pFlow::sphereParticlesKernels
 {
 
-using rpAcceleration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
+void addMassDiamInertiaProp(
+    deviceViewType1D<uint32>    shapeIndex,
+    deviceViewType1D<real>  	mass,
+	deviceViewType1D<real>  	diameter,
+	deviceViewType1D<real>  	I,
+    deviceViewType1D<uint32>  	propertyId,
+    pFlagTypeDevice 		    incld,
+    deviceViewType1D<real>  	src_mass,
+	deviceViewType1D<real>  	src_diameter,
+	deviceViewType1D<real>  	src_I,
+    deviceViewType1D<uint32>  	src_propertyId
+);
 
-template<typename IncludeFunctionType>
 void acceleration( 
 	realx3		g,
 	deviceViewType1D<real>  	mass,
 	deviceViewType1D<realx3>  	force,
 	deviceViewType1D<real>  	I,
 	deviceViewType1D<realx3>  	torque,
-	IncludeFunctionType 		incld,
+	pFlagTypeDevice 		    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]/mass[i] + g;
-				rAcc[i] = torque[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]/mass[i] + g;
-				rAcc[i] = torque[i]/I[i];
-			}
-		});
-
-	}
-	
-	Kokkos::fence();
-}
 
 }
 

src/Particles/SphereParticles/sphereShape/sphereShape.cpp

@@ -21,9 +21,9 @@ Licence:
 #include "sphereShape.hpp"
 
 
-bool pFlow::sphereShape::readDictionary()
+bool pFlow::sphereShape::readFromDictionary3()
 {
-	
+
 	diameters_ = getVal<realVector>("diameters");
 
 	if(diameters_.size() != numShapes() )
@@ -61,6 +61,11 @@ pFlow::sphereShape::sphereShape
 	shape(fileName, owner, prop)
 {
 
+	if(!readFromDictionary3())
+	{
+		fatalExit;
+		fatalErrorInFunction;
+	}
 }
 
 pFlow::real pFlow::sphereShape::maxBoundingSphere() const

src/Particles/SphereParticles/sphereShape/sphereShape.hpp

@@ -26,8 +26,6 @@ Licence:
 namespace pFlow
 {
 
-
-
 class sphereShape
 :
 	public shape
@@ -37,7 +35,7 @@ private:
 	// - diameter of spheres
 	realVector 			diameters_;
 	
-	bool readDictionary();
+	bool readFromDictionary3();
 
 protected: