From c202f9eaae316dcbd0f165b660e943cef9a7cf23 Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Mon, 20 Jan 2025 14:55:12 +0330 Subject: [PATCH 1/5] AB3, AB4 added, and AB2 modified --- .../AdamsBashforth2/AdamsBashforth2.hpp | 17 +- .../AdamsBashforth3/AdamsBashforth3.cpp | 166 ++++++++++++++--- .../AdamsBashforth3/AdamsBashforth3.hpp | 133 ++++++-------- .../AdamsBashforth4/AdamsBashforth4.cpp | 171 +++++++++++++++--- .../AdamsBashforth4/AdamsBashforth4.hpp | 121 ++++--------- src/Integration/CMakeLists.txt | 5 +- 6 files changed, 386 insertions(+), 227 deletions(-) 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 From 967bb753aa72270a3c0e028dfdf4b8d195aeb6dc Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Mon, 20 Jan 2025 15:37:48 +0330 Subject: [PATCH 2/5] adding non-linear contact force model for grains --- .../Models/contactForce/cGNonLinearCF.hpp | 307 ++++++++++++++++++ .../contactForce/cGRelativeLinearCF.hpp | 60 +--- .../Models/grainContactForceModels.hpp | 4 + .../grainInteractionsNonLinearModels.cpp | 6 + 4 files changed, 334 insertions(+), 43 deletions(-) create mode 100644 src/Interaction/Models/contactForce/cGNonLinearCF.hpp diff --git a/src/Interaction/Models/contactForce/cGNonLinearCF.hpp b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp new file mode 100644 index 00000000..a167f8fa --- /dev/null +++ b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp @@ -0,0 +1,307 @@ +/*------------------------------- 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 __cGNonLinearCF_hpp__ +#define __cGNonLinearCF_hpp__ + +#include "types.hpp" + +namespace pFlow::cfModels +{ + +template +class cGNonLinear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + struct cGNonLinearProperties + { + real Yeff_ = 1000000.0; + real Geff_ = 8000000.0; + real en_ = 1.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + cGNonLinearProperties(){} + + INLINE_FUNCTION_HD + cGNonLinearProperties(real Yeff, real Geff, real en, real mu ): + Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu) + {} + + INLINE_FUNCTION_HD + cGNonLinearProperties(const cGNonLinearProperties&)=default; + + INLINE_FUNCTION_HD + cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default; + + INLINE_FUNCTION_HD + ~cGNonLinearProperties() = default; + }; + +protected: + + using cGNonLinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + int32 addDissipationModel_; + + cGNonLinearArrayType nonlinearProperties_; + + bool readNonLinearDictionary(const dictionary& dict) + { + auto Yeff = dict.getVal("Yeff"); + auto Geff = dict.getVal("Geff"); + auto nu = dict.getVal("nu"); + auto en = dict.getVal("en"); + auto mu = dict.getVal("mu"); + + auto nElem = Yeff.size(); + + if(nElem != nu.size()) + { + fatalErrorInFunction<< + "sizes of Yeff("< prop("prop",nElem); + ForAll(i,Yeff) + { + prop[i] = {Yeff[i], Geff[i], en[i], mu[i]}; + } + + nonlinearProperties_.assign(prop); + + auto adm = dict.getVal("additionalDissipationModel"); + if(adm == "none") + { + addDissipationModel_ = 1; + } + else if(adm == "LU") + { + addDissipationModel_ = 2; + } + else if (adm == "GB") + { + addDissipationModel_ = 3; + } + else + { + addDissipationModel_ = 1; + } + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "cGNonLinearLimited"; + } + else + { + return "cGNonLinearNonLimited"; + } + return ""; + } + +public: + + TypeInfoNV(modelName()); + + INLINE_FUNCTION_HD + cGNonLinear(){} + + cGNonLinear( + int32 nMaterial, + const ViewType1D& rho, + const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + nonlinearProperties_("cGNonLinearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readNonLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + cGNonLinear(const cGNonLinear&) = default; + + INLINE_FUNCTION_HD + cGNonLinear(cGNonLinear&&) = default; + + INLINE_FUNCTION_HD + cGNonLinear& operator=(const cGNonLinear&) = default; + + INLINE_FUNCTION_HD + cGNonLinear& operator=(cGNonLinear&&) = default; + + + INLINE_FUNCTION_HD + ~cGNonLinear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, + const real Ri, + const real Rj, + const real cGFi, + const real cGFj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + const auto prop = nonlinearProperties_(propId_i,propId_j); + + const real f = 2.0/( 1.0/cGFi + 1.0/cGFj ); + + real vrn = dot(Vr, Nij); + realx3 Vt = Vr - vrn*Nij; + + history.overlap_t_ += Vt*dt; + + real mi = 3*Pi/4*pow(Ri,static_cast(3))*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,static_cast(3))*rho_[propId_j]; + real Reff = 1.0/(1/Ri + 1/Rj); + + real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff); + real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz); + + real en = prop.en_; + if (addDissipationModel_==2) + { + en = sqrt(1+((pow(prop.en_,2)-1)*f)); + } + else if (addDissipationModel_==3) + { + + en = exp((pow(f,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) )); + } + + real Kn = static_cast(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n); + + real ethan_ = -2.0*log(en)*sqrt(Kn)/ + sqrt(pow(log(en),2.0)+ pow(Pi,2.0)); + + FCn = ( - Kn*ovrlp_n - + sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast(0.25))*vrn)*Nij; + + FCt = (f*f)*(- static_cast(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_; + + real ft = length(FCt); + real ft_fric = prop.mu_ * length(FCn); + + // apply friction + if(ft > ft_fric) + { + if( length(history.overlap_t_) >0.0) + { + if constexpr (limited) + { + real kt = static_cast(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n); + FCt *= (ft_fric/ft); + history.overlap_t_ = - FCt/(f*f*kt); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + + } + else + { + FCt = 0.0; + } + } + + } + + +}; + +} //pFlow::CFModels + +#endif diff --git a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp index 48a77884..d542e70c 100755 --- a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp +++ b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp @@ -26,11 +26,6 @@ Licence: - - - - - namespace pFlow::cfModels { @@ -49,16 +44,15 @@ public: { real kn_ = 1000.0; real kt_ = 800.0; - real en_ = 0.0; - real ethat_ = 0.0; - real mu_ = 0.00001; + real en_ = 1.0; + real mu_ = 0.00001; INLINE_FUNCTION_HD linearProperties(){} INLINE_FUNCTION_HD - linearProperties(real kn, real kt, real en, real etha_t, real mu ): - kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) + linearProperties(real kn, real kt, real en, real mu ): + kn_(kn), kt_(kt), en_(en), mu_(mu) {} INLINE_FUNCTION_HD @@ -93,7 +87,6 @@ protected: auto kn = dict.getVal("kn"); auto kt = dict.getVal("kt"); auto en = dict.getVal("en"); - auto et = dict.getVal("et"); auto mu = dict.getVal("mu"); @@ -114,12 +107,6 @@ protected: return false; } - if(nElem != et.size()) - { - fatalErrorInFunction<< - "sizes of kn("< prop("prop", nElem); ForAll(i,kn) { - prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] }; + prop[i] = {kn[i], kt[i], en[i], mu[i] }; } linearProperties_.assign(prop); - auto adm = dict.getVal("additionalDissipationModel"); - + auto adm = dict.getVal("additionalDissipationModel"); if(adm == "none") { addDissipationModel_ = 1; @@ -287,22 +260,23 @@ public: real sqrt_meff = sqrt((mi*mj)/(mi+mj)); - - if (addDissipationModel_==2) + real en = prop.en_; + if (addDissipationModel_==2) { - prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_)); + en = sqrt(1+((pow(prop.en_,2)-1)*f_)); } - else if (addDissipationModel_==3) + else if (addDissipationModel_==3) { - auto pie =3.14; - prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); + + en = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) )); } - real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/ - sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); + + real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/ + sqrt(pow(log(en),2.0)+ pow(Pi,2.0)); - FCn = ( -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; - FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) * prop.ethat_*Vt); + FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; + FCt = ( -f_*prop.kt_ * history.overlap_t_); diff --git a/src/Interaction/Models/grainContactForceModels.hpp b/src/Interaction/Models/grainContactForceModels.hpp index 908117f5..1b4d4cbe 100755 --- a/src/Interaction/Models/grainContactForceModels.hpp +++ b/src/Interaction/Models/grainContactForceModels.hpp @@ -23,6 +23,7 @@ Licence: #include "cGAbsoluteLinearCF.hpp" #include "cGRelativeLinearCF.hpp" +#include "cGNonLinearCF.hpp" #include "grainRolling.hpp" @@ -38,6 +39,9 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling>; using nonLimitedCGRelativeLinearGrainRolling = grainRolling>; +using limitedCGNonLinearGrainRolling = grainRolling>; +using nonLimitedCGNonLinearGrainRolling = grainRolling>; + } diff --git a/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp index 462b5eb3..42f0b469 100755 --- a/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp +++ b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp @@ -41,3 +41,9 @@ Licence: createBoundaryGrainInteraction(ForceModel, GeomModel) +createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry); + +createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry); + +createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry); +createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry); \ No newline at end of file From a32786eb8a6268feed3d8d0bf476574813bcc69a Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Mon, 20 Jan 2025 15:39:17 +0330 Subject: [PATCH 3/5] particlePhasicFlow bug-fix when flag --set-fields-only is used and no shpaeName is set --- .../particlesPhasicFlow.cpp | 50 ++++++++++--------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp index 637fcf65..9729e27a 100755 --- a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp +++ b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp @@ -181,32 +181,36 @@ int main( int argc, char* argv[] ) &Control.caseSetup() ); - auto& shapeName = Control.time().template lookupObject("shapeName"); - - REPORT(0)<< "Converting shapeName field to shapeIndex field"<("shapeName"); + + REPORT(0)<< "Converting shapeName field to shapeIndex field"< Date: Mon, 20 Jan 2025 15:43:56 +0330 Subject: [PATCH 4/5] bug fixes for build with float in cuda --- CMakeLists.txt | 2 +- src/phasicFlow/dictionary/dictionary.hpp | 2 +- src/phasicFlow/repository/Time/timeControl.cpp | 2 +- src/phasicFlow/structuredData/domain/simulationDomain.cpp | 2 +- src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp | 2 +- src/phasicFlow/types/basicTypes/bTypesFunctions.cpp | 1 - src/phasicFlow/types/basicTypes/bTypesFunctions.hpp | 2 +- 7 files changed, 6 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fbb2eaf2..5bdcde81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,7 @@ project(phasicFlow VERSION 1.0 ) set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE) set(CMAKE_CXX_STANDARD_REQUIRED True) set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE) -set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE) +set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE) set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE) mark_as_advanced(FORCE var BUILD_SHARED_LIBS) diff --git a/src/phasicFlow/dictionary/dictionary.hpp b/src/phasicFlow/dictionary/dictionary.hpp index e64dfec1..0a7da6d5 100644 --- a/src/phasicFlow/dictionary/dictionary.hpp +++ b/src/phasicFlow/dictionary/dictionary.hpp @@ -431,7 +431,7 @@ T dictionary::getValOrSet template T dictionary::getValOrSetMax(const word& keyword, const T& setMaxVal)const { - return max(getValOrSet(keyword, setMaxVal), setMaxVal); + return max(getValOrSet(keyword, setMaxVal), static_cast(setMaxVal)); } template diff --git a/src/phasicFlow/repository/Time/timeControl.cpp b/src/phasicFlow/repository/Time/timeControl.cpp index b6305545..6ccc4086 100644 --- a/src/phasicFlow/repository/Time/timeControl.cpp +++ b/src/phasicFlow/repository/Time/timeControl.cpp @@ -167,7 +167,7 @@ void pFlow::timeControl::checkForOutputToFile() lastSaved_ = currentTime_; save = true; } - else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) ) + else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), static_cast(0.5 *dt_)) ) { lastSaved_ = currentTime_; save = true; diff --git a/src/phasicFlow/structuredData/domain/simulationDomain.cpp b/src/phasicFlow/structuredData/domain/simulationDomain.cpp index eaf55965..7621e79f 100644 --- a/src/phasicFlow/structuredData/domain/simulationDomain.cpp +++ b/src/phasicFlow/structuredData/domain/simulationDomain.cpp @@ -31,7 +31,7 @@ bool pFlow::simulationDomain::prepareBoundaryDicts() real neighborLength = boundaries.getVal("neighborLength"); real boundaryExtntionLengthRatio = - boundaries.getValOrSetMax("boundaryExtntionLengthRatio", 0.1); + boundaries.getValOrSetMax("boundaryExtntionLengthRatio", static_cast(0.1)); uint32 updateInterval = boundaries.getValOrSetMax("updateInterval", 1u); diff --git a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp index d9f41a38..89e497fb 100644 --- a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp +++ b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp @@ -28,7 +28,7 @@ pFlow::infinitePlane::infinitePlane ) { auto ln = cross(p2-p1, p3-p1); - + if( equal(ln.length(),0.0) ) { fatalErrorInFunction<< diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp b/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp index c99eb4bd..dedf002d 100755 --- a/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp +++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp @@ -21,7 +21,6 @@ Licence: #include #include #include - #include "bTypesFunctions.hpp" pFlow::int32 diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp b/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp index 6efe0e28..5ad36fbb 100755 --- a/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp +++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp @@ -26,7 +26,7 @@ Licence: #define __bTypesFunctions_hpp__ #include "builtinTypes.hpp" -//#include "math.hpp" +#include "math.hpp" #include "numericConstants.hpp" #include "pFlowMacros.hpp" From 2ec3dfdc7ed30a901a1a73f23336bf2dddcb3346 Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Mon, 20 Jan 2025 21:02:50 +0330 Subject: [PATCH 5/5] global damping is added to the code --- .../AdamsBashforth2/AdamsBashforth2.cpp | 17 +++-- .../AdamsBashforth2/AdamsBashforth2.hpp | 3 +- .../AdamsBashforth3/AdamsBashforth3.cpp | 21 +++--- .../AdamsBashforth3/AdamsBashforth3.hpp | 3 +- .../AdamsBashforth4/AdamsBashforth4.cpp | 21 +++--- .../AdamsBashforth4/AdamsBashforth4.hpp | 3 +- src/Integration/integration/integration.hpp | 2 +- src/Particles/CMakeLists.txt | 3 + .../dynamicPointStructure.cpp | 17 ++++- .../dynamicPointStructure.hpp | 6 ++ src/Particles/globalDamping/globalDamping.cpp | 74 +++++++++++++++++++ src/Particles/globalDamping/globalDamping.hpp | 65 ++++++++++++++++ .../repository/Time/baseTimeControl.cpp | 44 +++++++++++ .../repository/Time/baseTimeControl.hpp | 6 ++ 14 files changed, 255 insertions(+), 30 deletions(-) create mode 100644 src/Particles/globalDamping/globalDamping.cpp create mode 100644 src/Particles/globalDamping/globalDamping.hpp diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp b/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp index 7969b184..09733eb4 100644 --- a/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp @@ -37,7 +37,8 @@ bool intAllActive( real dt, realx3Field_D& y, realx3PointField_D& dy, - realx3PointField_D& dy1) + realx3PointField_D& dy1, + real damping = 1.0) { auto d_dy = dy.deviceView(); @@ -49,7 +50,7 @@ bool intAllActive( "AdamsBashforth2::correct", rpIntegration (activeRng.start(), activeRng.end()), LAMBDA_HD(uint32 i){ - d_y[i] += dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i]); + d_y[i] = damping*(d_dy[i] + dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i])); d_dy1[i] = d_dy[i]; }); Kokkos::fence(); @@ -62,7 +63,8 @@ bool intScattered real dt, realx3Field_D& y, realx3PointField_D& dy, - realx3PointField_D& dy1 + realx3PointField_D& dy1, + real damping = 1.0 ) { @@ -78,7 +80,7 @@ bool intScattered LAMBDA_HD(uint32 i){ if( activeP(i)) { - d_y[i] += dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i]); + d_y[i] = damping*(d_y[i] + dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i])); d_dy1[i] = d_dy[i]; } }); @@ -142,18 +144,19 @@ bool pFlow::AdamsBashforth2::correct ( real dt, realx3PointField_D& y, - realx3PointField_D& dy + realx3PointField_D& dy, + real damping ) { auto& dy1l = dy1(); bool success = false; if(dy1l.isAllActive()) { - success = intAllActive(dt, y.field(), dy, dy1l); + success = intAllActive(dt, y.field(), dy, dy1l, damping); } else { - success = intScattered(dt, y.field(), dy, dy1l); + success = intScattered(dt, y.field(), dy, dy1l, damping); } success = success && boundaryList_.correct(dt, y, dy); diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp index e0e92413..91b5902b 100644 --- a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp @@ -109,7 +109,8 @@ public: bool correct( real dt, realx3PointField_D& y, - realx3PointField_D& dy) override; + realx3PointField_D& dy, + real damping = 1.0) override; bool correctPStruct( real dt, diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp index e7cdf0d0..e8198512 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp @@ -35,7 +35,8 @@ bool intAllActive( realx3Field_D& y, realx3PointField_D& dy, realx3PointField_D& dy1, - realx3PointField_D& dy2) + realx3PointField_D& dy2, + real damping = 1.0) { auto d_dy = dy.deviceView(); @@ -48,9 +49,9 @@ bool intAllActive( "AdamsBashforth3::correct", rpIntegration (activeRng.start(), activeRng.end()), LAMBDA_HD(uint32 i){ - d_y[i] += dt*( static_cast(23.0 / 12.0) * d_dy[i] + d_y[i] = damping*( 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]); + + static_cast(5.0 / 12.0) * d_dy2[i]) ); d_dy2[i] = d_dy1[i]; d_dy1[i] = d_dy[i]; @@ -67,7 +68,8 @@ bool intScattered realx3Field_D& y, realx3PointField_D& dy, realx3PointField_D& dy1, - realx3PointField_D& dy2 + realx3PointField_D& dy2, + real damping = 1.0 ) { @@ -84,9 +86,9 @@ bool intScattered LAMBDA_HD(uint32 i){ if( activeP(i)) { - d_y[i] += dt*( static_cast(23.0 / 12.0) * d_dy[i] + d_y[i] = damping * (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]); + + static_cast(5.0 / 12.0) * d_dy2[i])); d_dy2[i] = d_dy1[i]; d_dy1[i] = d_dy[i]; @@ -141,18 +143,19 @@ bool pFlow::AdamsBashforth3::correct ( real dt, realx3PointField_D& y, - realx3PointField_D& dy + realx3PointField_D& dy, + real damping ) { bool success = false; if(y.isAllActive()) { - success = intAllActive(dt, y.field(), dy, dy1(), dy2()); + success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping); } else { - success = intScattered(dt, y.field(), dy, dy1(), dy2()); + success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping); } success = success && boundaryList().correct(dt, y, dy); diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp index fa21543f..e9f71c69 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp @@ -99,7 +99,8 @@ public: bool correct( real dt, realx3PointField_D& y, - realx3PointField_D& dy) override; + realx3PointField_D& dy, + real damping = 1.0) override; bool correctPStruct( real dt, diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp index 5764a20d..b6cc481b 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp +++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp @@ -36,7 +36,8 @@ bool intAllActive( realx3PointField_D& dy, realx3PointField_D& dy1, realx3PointField_D& dy2, - realx3PointField_D& dy3) + realx3PointField_D& dy3, + real damping = 1.0) { auto d_dy = dy.deviceView(); @@ -50,12 +51,12 @@ bool intAllActive( "AdamsBashforth3::correct", rpIntegration (activeRng.start(), activeRng.end()), LAMBDA_HD(uint32 i){ - d_y[i] += dt*( + d_y[i] = damping *(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]; @@ -72,7 +73,8 @@ bool intScattered realx3PointField_D& dy, realx3PointField_D& dy1, realx3PointField_D& dy2, - realx3PointField_D& dy3 + realx3PointField_D& dy3, + real damping = 1.0 ) { @@ -90,12 +92,12 @@ bool intScattered LAMBDA_HD(uint32 i){ if( activeP(i)) { - d_y[i] += dt*( + d_y[i] = damping* ( 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]; @@ -147,18 +149,19 @@ bool pFlow::AdamsBashforth4::correct ( real dt, realx3PointField_D& y, - realx3PointField_D& dy + realx3PointField_D& dy, + real damping ) { bool success = false; if(y.isAllActive()) { - success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3()); + success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping); } else { - success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3()); + success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping); } success = success && boundaryList().correct(dt, y, dy); diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp index 02e91611..32d82c22 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp +++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp @@ -96,7 +96,8 @@ public: bool correct( real dt, realx3PointField_D& y, - realx3PointField_D& dy) override; + realx3PointField_D& dy, + real damping = 1.0) override; bool correctPStruct( real dt, diff --git a/src/Integration/integration/integration.hpp b/src/Integration/integration/integration.hpp index 2c95813e..94f91f70 100644 --- a/src/Integration/integration/integration.hpp +++ b/src/Integration/integration/integration.hpp @@ -146,7 +146,7 @@ public: /// Correction/main integration step virtual - bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0; + bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0; virtual bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0; diff --git a/src/Particles/CMakeLists.txt b/src/Particles/CMakeLists.txt index 7cc91295..b0728d5a 100644 --- a/src/Particles/CMakeLists.txt +++ b/src/Particles/CMakeLists.txt @@ -6,6 +6,9 @@ particles/shape/shape.cpp particles/particles.cpp particles/particleIdHandler/particleIdHandler.cpp particles/regularParticleIdHandler/regularParticleIdHandler.cpp + +globalDamping/globalDamping.cpp + SphereParticles/sphereShape/sphereShape.cpp SphereParticles/sphereParticles/sphereParticles.cpp SphereParticles/sphereParticles/sphereParticlesKernels.cpp diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp index f0875341..21289928 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp @@ -21,6 +21,7 @@ Licence: #include "dynamicPointStructure.hpp" #include "systemControl.hpp" + pFlow::dynamicPointStructure::dynamicPointStructure ( systemControl& control @@ -77,6 +78,9 @@ pFlow::dynamicPointStructure::dynamicPointStructure fatalExit; } + REPORT(1)<<"Reading globalDamping dictionary ..."<(control); + } @@ -101,6 +105,16 @@ bool pFlow::dynamicPointStructure::iterate() return correct(dt, acc);*/ } +bool pFlow::dynamicPointStructure::afterIteration() +{ + //const auto ti = TimeInfo(); + + auto succs = pointStructure::afterIteration(); + //velDamping_().applyDamping(ti, velocity_); + + return succs; +} + bool pFlow::dynamicPointStructure::predict( real dt, realx3PointField_D &acceleration) @@ -119,10 +133,11 @@ bool pFlow::dynamicPointStructure::correct ) { //auto& pos = pStruct().pointPosition(); + const auto ti = TimeInfo(); if(!integrationPos_().correctPStruct(dt, *this, velocity_) )return false; - if(!integrationVel_().correct(dt, velocity_, acceleration))return false; + if(!integrationVel_().correct(dt, velocity_, acceleration, velDamping_().dampingFactor(ti)))return false; return true; } diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp index f82f1dca..81079a66 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp @@ -26,11 +26,13 @@ Licence: #include "pointFields.hpp" #include "integration.hpp" #include "uniquePtr.hpp" +#include "globalDamping.hpp" namespace pFlow { class systemControl; +//class globalDamping; class dynamicPointStructure : @@ -44,6 +46,8 @@ private: uniquePtr integrationVel_ = nullptr; + uniquePtr velDamping_ = nullptr; + Timer velocityUpdateTimer_; /// @brief integration method for velocity and position @@ -88,6 +92,8 @@ public: /// @brief This is called in time loop. Perform the main calculations /// when the component should evolve along time. bool iterate() override; + + bool afterIteration()override; /// prediction step (if any), is called in beforeIteration bool predict(real dt, realx3PointField_D& acceleration); diff --git a/src/Particles/globalDamping/globalDamping.cpp b/src/Particles/globalDamping/globalDamping.cpp new file mode 100644 index 00000000..633d70fb --- /dev/null +++ b/src/Particles/globalDamping/globalDamping.cpp @@ -0,0 +1,74 @@ +/*------------------------------- 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 "globalDamping.hpp" + + +pFlow::globalDamping::globalDamping(const systemControl& control) +: + timeControl_(control.settingsDict().subDict("globalDamping"), control.time().dt(), "damping") +{ + const dictionary& dict = control.settingsDict().subDict("globalDamping"); + + dampingFactor_ = dict.getValOrSetMin("dampingFactor", static_cast(1.0)); + + dampingFactor_ = max( dampingFactor_ , static_cast(0.01)); + + performDamping_ = !equal(dampingFactor_, static_cast(1.0)); + + if( performDamping_ ) + REPORT(2)<<"Global damping "<("timeControl"); + if(tControl == "timeStep") + isTimeStep_ = true; + else if( tControl == "runTime" || + tControl == "simulationTime") + isTimeStep_ = false; + else + { + fatalErrorInFunction<< + "bad input for intervalControl in dictionary "<< dict.globalName()<("startTime", defStartTime)); + auto endTime = (dict.getValOrSet("endTime", largeValue)); + auto interval = dict.getValOrSet(intervalWord, defInterval); + rRange_ = realStridedRange(startTime, endTime, interval); + + } + else + { + auto startTime = (dict.getValOrSet("startTime", 0)); + auto endTime = (dict.getValOrSet("endTime", largestPosInt32)); + auto interval = dict.getValOrSet(intervalWord, 1); + iRange_ = int32StridedRagne(startTime, endTime, interval); + } +} + pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix) : isTimeStep_(true), diff --git a/src/phasicFlow/repository/Time/baseTimeControl.hpp b/src/phasicFlow/repository/Time/baseTimeControl.hpp index 5a7634d5..cb0e5d3e 100644 --- a/src/phasicFlow/repository/Time/baseTimeControl.hpp +++ b/src/phasicFlow/repository/Time/baseTimeControl.hpp @@ -46,6 +46,12 @@ public: real defStartTime = 0.0 ); + baseTimeControl( + const dictionary& dict, + const real defInterval, + const word& intervalPrefix="", + const real defStartTime=0.0); + baseTimeControl( int32 start, int32 end,