AB3, AB4 added, and AB2 modified

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -45,6 +45,8 @@ private:
 
 	friend class processorAB2BoundaryIntegration;
 
+protected:
+	
 	const auto& dy1()const
 	{
 		return static_cast<const realx3PointField_D&>(*this);
@@ -55,6 +57,11 @@ private:
 		return static_cast<realx3PointField_D&>(*this);
 	}
 
+	boundaryIntegrationList& boundaryList()
+	{
+		return boundaryList_;
+	}
+	
 public:
 
 	/// Class info
@@ -70,7 +77,7 @@ public:
 			const realx3Field_D& initialValField);
 
 		/// Destructor 
-		~AdamsBashforth2()final = default;
+		~AdamsBashforth2()override = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -102,12 +109,12 @@ public:
 		bool correct(
 			real dt, 
 			realx3PointField_D& y, 
-			realx3PointField_D& dy) final;
+			realx3PointField_D& dy) override;
 
 		bool correctPStruct(
 			real dt, 
 			pointStructure& pStruct, 
-			realx3PointField_D& vel) final;
+			realx3PointField_D& vel) override;
 			
 		
 		/*bool hearChanges
@@ -121,9 +128,9 @@ public:
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
-			const realx3Vector& y) final;
+			const realx3Vector& y) override;
 
-		bool needSetInitialVals()const final
+		bool needSetInitialVals()const override
 		{
 			return false;
 		}

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -20,62 +20,168 @@ Licence:
 
 #include "AdamsBashforth3.hpp"
 
+namespace pFlow
+{
+
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
+
+bool intAllActive(
+	real dt, 
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2)
+{
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto d_dy2= dy2.deviceView();
+	auto activeRng = dy1.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+			
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
+
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth2::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+			
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
 //const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
 
 pFlow::AdamsBashforth3::AdamsBashforth3
 (
 	const word& baseName,
-	repository& owner,
+	pointStructure& pStruct,
-	const pointStructure& pStruct,
+	const word& method,
-	const word& method
+	const realx3Field_D& initialValField
 )
 :
-	integration(baseName, owner, pStruct, method),
+	AdamsBashforth2(baseName, pStruct, method, initialValField),
-	history_(
+	dy2_
-		owner.emplaceObject<pointField<VectorSingle,AB3History>>(
+	(
-			objectFile(
+		objectFile
-				groupNames(baseName,"AB3History"),
+		(
-				"",
+			groupNames(baseName,"dy2"),
-				objectFile::READ_IF_PRESENT,
+			pStruct.time().integrationFolder(),
-				objectFile::WRITE_ALWAYS),
+			objectFile::READ_IF_PRESENT,
-			pStruct,
+			objectFile::WRITE_ALWAYS
-			AB3History({zero3,zero3})))
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
 
 {
 	
 }
 
-bool pFlow::AdamsBashforth3::predict
+void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
 {
-
+	AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
-	return true;
+	dy2_.updateBoundariesSlaveToMasterIfRequested();
 }
 
+
+
 bool pFlow::AdamsBashforth3::correct
 (
 	real dt, 
-	realx3Vector_D& y,
+	realx3PointField_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy
 )
 {
 	
-	if(this->pStruct().allActive())
+	bool success = false;
+	if(y.isAllActive())
 	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2());
 	}
 	else
 	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+		success = intScattered(dt, y.field(), dy, dy1(), dy2());
 	}
 
-	return true;
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
 }
 
+bool pFlow::AdamsBashforth3::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+
 bool pFlow::AdamsBashforth3::setInitialVals(
 	const int32IndexContainer& newIndices,
 	const realx3Vector& y)
@@ -83,7 +189,7 @@ bool pFlow::AdamsBashforth3::setInitialVals(
 	return true;
 }
 
-bool pFlow::AdamsBashforth3::intAll(
+/*bool pFlow::AdamsBashforth3::intAll(
 	real dt, 
 	realx3Vector_D& y, 
 	realx3Vector_D& dy, 
@@ -106,4 +212,4 @@ bool pFlow::AdamsBashforth3::intAll(
 	Kokkos::fence();
 
 	return true;	
-}
+}*/

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -22,48 +22,14 @@ Licence:
 #define __AdamsBashforth3_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth2.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB3History
-{
-	TypeInfoNV("AB3History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB3History& ab3)
-{
-	str.readBegin("AB3History");
-
-    str >> ab3.dy1_;
-    str >> ab3.dy2_;
-    
-    str.readEnd("AB3History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB3History& ab3)
-{
-	str << token::BEGIN_LIST << ab3.dy1_
-	    << token::SPACE << ab3.dy2_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Third order Adams-Bashforth integration method for solving ODE
@@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
  */
 class AdamsBashforth3
 :
-	public integration
+	public AdamsBashforth2
 {
+private:
+	
+	friend class processorAB3BoundaryIntegration;
+
+	realx3PointField_D 		dy2_;
+
 protected:
 
-	/// Integration history
+	const auto& dy2()const
-	pointField<VectorSingle,AB3History>& history_;
+	{
+		return dy2_;
+	}
+
+	auto& dy2()
+	{
+		return dy2_;
+	}
 	
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -96,17 +69,13 @@ public:
 		/// Construct from components
 		AdamsBashforth3(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth3>(*this);
-		}
 
 		/// Destructor 
-		virtual ~AdamsBashforth3()=default;
+		~AdamsBashforth3() override =default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -117,14 +86,35 @@ public:
 
 	// - Methods
 
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, 
+		/// return integration method 
-			realx3Vector_D & y,
+		word method()const override
-			realx3Vector_D& dy) override;
+		{
+			return "AdamsBashforth3";
+		}
+
+		
+
+		bool correct(
+			real dt, 
+			realx3PointField_D& y, 
+			realx3PointField_D& dy) override;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
+			
+		
+		/*bool hearChanges
+		(
+			real t,
+			real dt,
+			uint32 iter,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -135,25 +125,10 @@ public:
 			return false;
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
-
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
-
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth3::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -181,7 +156,7 @@ bool pFlow::AdamsBashforth3::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -20,60 +20,171 @@ Licence:
 
 #include "AdamsBashforth4.hpp"
 
+namespace pFlow
+{
 
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
+
+bool intAllActive(
+	real dt, 
+	realx3Field_D& y, 
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3)
+{
+
+	auto d_dy = dy.deviceView();
+	auto d_y  = y.deviceView();
+	auto d_dy1= dy1.deviceView();
+	auto d_dy2= dy2.deviceView();
+	auto d_dy3= dy3.deviceView();
+	auto activeRng = dy1.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_y[i] += dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
+				);
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = d_dy1[i];
+			d_dy1[i] = d_dy[i];
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+bool intScattered
+(
+	real dt, 
+	realx3Field_D& y,
+	realx3PointField_D& dy,
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth2::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
+				);
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
 
 pFlow::AdamsBashforth4::AdamsBashforth4
 (
 	const word& baseName,
-	repository& owner,
+	pointStructure& pStruct,
-	const pointStructure& pStruct,
+	const word& method,
-	const word& method
+	const realx3Field_D& initialValField
 )
 :
-	integration(baseName, owner, pStruct, method),
+	AdamsBashforth3(baseName, pStruct, method, initialValField),
-	history_(
+	dy3_
-		owner.emplaceObject<pointField<VectorSingle,AB4History>>(
+	(
-			objectFile(
+		objectFile
-				groupNames(baseName,"AB4History"),
+		(
-				"",
+			groupNames(baseName,"dy3"),
-				objectFile::READ_IF_PRESENT,
+			pStruct.time().integrationFolder(),
-				objectFile::WRITE_ALWAYS),
+			objectFile::READ_IF_PRESENT,
-			pStruct,
+			objectFile::WRITE_ALWAYS
-			AB4History({zero3,zero3, zero3})))
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+	
 
 {
 
 }
 
-bool pFlow::AdamsBashforth4::predict
+void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
 {
-
+	AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
-	return true;
+	dy3_.updateBoundariesSlaveToMasterIfRequested();
 }
 
 bool pFlow::AdamsBashforth4::correct
 (
 	real dt, 
-	realx3Vector_D& y,
+	realx3PointField_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy
 )
 {
 	
-	if(this->pStruct().allActive())
+	bool success = false;
+	if(y.isAllActive())
 	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3());
 	}
 	else
 	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3());
 	}
 
-	return true;
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth4::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
 }
 
 bool pFlow::AdamsBashforth4::setInitialVals(
@@ -83,7 +194,7 @@ bool pFlow::AdamsBashforth4::setInitialVals(
 	return true;
 }
 
-bool pFlow::AdamsBashforth4::intAll(
+/*bool pFlow::AdamsBashforth4::intAll(
 	real dt, 
 	realx3Vector_D& y, 
 	realx3Vector_D& dy, 
@@ -112,4 +223,4 @@ bool pFlow::AdamsBashforth4::intAll(
 	Kokkos::fence();
 
 	return true;	
-}
+}*/

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -22,53 +22,14 @@ Licence:
 #define __AdamsBashforth4_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth3.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB4History
-{
-	TypeInfoNV("AB4History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB4History& ab4)
-{
-	str.readBegin("AB4History");
-
-    str >> ab4.dy1_;
-    str >> ab4.dy2_;
-    str >> ab4.dy3_;
-    
-    str.readEnd("AB4History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB4History& ab4)
-{
-	str << token::BEGIN_LIST << ab4.dy1_
-	    << token::SPACE << ab4.dy2_
-	    << token::SPACE << ab4.dy3_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
-
 /**
  * Fourth order Adams-Bashforth integration method for solving ODE
  * 
@@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
  */
 class AdamsBashforth4
 :
-	public integration
+	public AdamsBashforth3
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy3_;
+
 protected:
 
-	/// Integration history
+	const auto& dy3()const
-	pointField<VectorSingle,AB4History>& history_;
+	{
+		return dy3_;
+	}
 
-	/// Range policy for integration kernel
+	auto& dy3()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return dy3_;
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -100,17 +67,14 @@ public:
 		/// Construct from components
 		AdamsBashforth4(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField);
+		
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth4>(*this);
-		}
 
 		/// Destructor
-		virtual ~AdamsBashforth4()=default;
+		~AdamsBashforth4() override =default;
 
 		/// Add a this to the virtual constructor table 
 		add_vCtor(
@@ -121,15 +85,23 @@ public:
 
 	// - Methods
 		
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
+
-			realx3Vector_D & UNUSED(y),
+		/// return integration method 
-			realx3Vector_D& UNUSED(dy)) override;
+		word method()const override
+		{
+			return "AdamsBashforth4";
+		}
 
 		bool correct(
 			real dt, 
-			realx3Vector_D & y,
+			realx3PointField_D& y, 
-			realx3Vector_D& dy) override;
+			realx3PointField_D& dy) override;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -140,25 +112,12 @@ public:
 			return false;
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
 		
-		/// Integrate on active points in the active range 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth4::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -191,7 +150,7 @@ bool pFlow::AdamsBashforth4::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/CMakeLists.txt

@@ -4,9 +4,10 @@ integration/integration.cpp
 boundaries/boundaryIntegration.cpp
 boundaries/boundaryIntegrationList.cpp 
 AdamsBashforth2/AdamsBashforth2.cpp
+AdamsBashforth3/AdamsBashforth3.cpp
+AdamsBashforth4/AdamsBashforth4.cpp
 #AdamsBashforth5/AdamsBashforth5.cpp
-#AdamsBashforth4/AdamsBashforth4.cpp
+
-#AdamsBashforth3/AdamsBashforth3.cpp
 #AdamsMoulton3/AdamsMoulton3.cpp
 #AdamsMoulton4/AdamsMoulton4.cpp
 #AdamsMoulton5/AdamsMoulton5.cpp