diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp index e8abc026..e0e92413 100644 --- a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp @@ -45,6 +45,8 @@ private: friend class processorAB2BoundaryIntegration; +protected: + const auto& dy1()const { return static_cast(*this); @@ -54,6 +56,11 @@ private: { return static_cast(*this); } + + boundaryIntegrationList& boundaryList() + { + return boundaryList_; + } public: @@ -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; } diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp index f8ec6b5a..e7cdf0d0 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp +++ b/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::IndexType + >; + +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(23.0 / 12.0) * d_dy[i] + - static_cast(16.0 / 12.0) * d_dy1[i] + + static_cast(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(23.0 / 12.0) * d_dy[i] + - static_cast(16.0 / 12.0) * d_dy1[i] + + static_cast(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, - const pointStructure& pStruct, - const word& method + pointStructure& pStruct, + const word& method, + const realx3Field_D& initialValField ) : - integration(baseName, owner, pStruct, method), - history_( - owner.emplaceObject>( - objectFile( - groupNames(baseName,"AB3History"), - "", - objectFile::READ_IF_PRESENT, - objectFile::WRITE_ALWAYS), - pStruct, - AB3History({zero3,zero3}))) + AdamsBashforth2(baseName, pStruct, method, initialValField), + dy2_ + ( + objectFile + ( + groupNames(baseName,"dy2"), + pStruct.time().integrationFolder(), + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct, + zero3, + zero3 + ) { } -bool pFlow::AdamsBashforth3::predict -( - real UNUSED(dt), - realx3Vector_D& UNUSED(y), - realx3Vector_D& UNUSED(dy) -) +void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested() { - - return true; + AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested(); + dy2_.updateBoundariesSlaveToMasterIfRequested(); } + + bool pFlow::AdamsBashforth3::correct ( - real dt, - realx3Vector_D& y, - realx3Vector_D& dy + real dt, + realx3PointField_D& y, + 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; -} \ No newline at end of file +}*/ \ No newline at end of file diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp index 279fdcda..fa21543f 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp +++ b/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 { -protected: +private: - /// Integration history - pointField& history_; + friend class processorAB3BoundaryIntegration; - /// Range policy for integration kernel - using rpIntegration = Kokkos::RangePolicy< - DefaultExecutionSpace, - Kokkos::Schedule, - Kokkos::IndexType - >; + realx3PointField_D dy2_; + +protected: + + const auto& dy2()const + { + return dy2_; + } + + auto& dy2() + { + return dy2_; + } + public: @@ -96,17 +69,13 @@ public: /// Construct from components AdamsBashforth3( const word& baseName, - repository& owner, - const pointStructure& pStruct, - const word& method); + pointStructure& pStruct, + const word& method, + const realx3Field_D& initialValField); - uniquePtr clone()const override - { - return makeUnique(*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( - real UNUSED(dt), - realx3Vector_D & UNUSED(y), - realx3Vector_D& UNUSED(dy)) override; + void updateBoundariesSlaveToMasterIfRequested()override; - bool correct(real dt, - realx3Vector_D & y, - realx3Vector_D& dy) override; + /// return integration method + word method()const 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 - bool intRange( - real dt, - realx3Vector_D& y, - realx3Vector_D& dy, - activeFunctor activeP ); - }; -template +/*template bool pFlow::AdamsBashforth3::intRange( real dt, realx3Vector_D& y, @@ -181,7 +156,7 @@ bool pFlow::AdamsBashforth3::intRange( Kokkos::fence(); return true; -} +}*/ } // pFlow diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp index 23bad29b..5764a20d 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp +++ b/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::IndexType + >; + +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(55.0 / 24.0) * d_dy[i] + - static_cast(59.0 / 24.0) * d_dy1[i] + + static_cast(37.0 / 24.0) * d_dy2[i] + - static_cast( 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(55.0 / 24.0) * d_dy[i] + - static_cast(59.0 / 24.0) * d_dy1[i] + + static_cast(37.0 / 24.0) * d_dy2[i] + - static_cast( 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, - const pointStructure& pStruct, - const word& method + pointStructure& pStruct, + const word& method, + const realx3Field_D& initialValField ) : - integration(baseName, owner, pStruct, method), - history_( - owner.emplaceObject>( - objectFile( - groupNames(baseName,"AB4History"), - "", - objectFile::READ_IF_PRESENT, - objectFile::WRITE_ALWAYS), - pStruct, - AB4History({zero3,zero3, zero3}))) + AdamsBashforth3(baseName, pStruct, method, initialValField), + dy3_ + ( + objectFile + ( + groupNames(baseName,"dy3"), + pStruct.time().integrationFolder(), + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + pStruct, + zero3, + zero3 + ) + { } -bool pFlow::AdamsBashforth4::predict -( - real UNUSED(dt), - realx3Vector_D& UNUSED(y), - realx3Vector_D& UNUSED(dy) -) +void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested() { - - return true; + AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested(); + dy3_.updateBoundariesSlaveToMasterIfRequested(); } bool pFlow::AdamsBashforth4::correct ( - real dt, - realx3Vector_D& y, - realx3Vector_D& dy + real dt, + realx3PointField_D& y, + 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; -} \ No newline at end of file +}*/ \ No newline at end of file diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp index cf1f50d8..02e91611 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp +++ b/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 - pointField& history_; + const auto& dy3()const + { + return dy3_; + } - /// Range policy for integration kernel - using rpIntegration = Kokkos::RangePolicy< - DefaultExecutionSpace, - Kokkos::Schedule, - Kokkos::IndexType - >; + auto& dy3() + { + return dy3_; + } public: @@ -100,17 +67,14 @@ public: /// Construct from components AdamsBashforth4( const word& baseName, - repository& owner, - const pointStructure& pStruct, - const word& method); + pointStructure& pStruct, + const word& method, + const realx3Field_D& initialValField); + - uniquePtr clone()const override - { - return makeUnique(*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( - real UNUSED(dt), - realx3Vector_D & UNUSED(y), - realx3Vector_D& UNUSED(dy)) override; + void updateBoundariesSlaveToMasterIfRequested()override; + + /// return integration method + word method()const override + { + return "AdamsBashforth4"; + } bool correct( real dt, - realx3Vector_D & y, - realx3Vector_D& dy) override; + realx3PointField_D& y, + 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 - bool intRange( - real dt, - realx3Vector_D& y, - realx3Vector_D& dy, - activeFunctor activeP ); + }; -template +/*template bool pFlow::AdamsBashforth4::intRange( real dt, realx3Vector_D& y, @@ -191,7 +150,7 @@ bool pFlow::AdamsBashforth4::intRange( Kokkos::fence(); return true; -} +}*/ } // pFlow diff --git a/src/Integration/CMakeLists.txt b/src/Integration/CMakeLists.txt index 32c3c558..ee8dfcfc 100644 --- a/src/Integration/CMakeLists.txt +++ b/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