diff --git a/CMakeLists.txt b/CMakeLists.txt index 2c45ef0b..f9740a90 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,7 +77,7 @@ add_subdirectory(solvers) add_subdirectory(utilities) add_subdirectory(DEMSystems) -#add_subdirectory(test_newFeatures) +#add_subdirectory(testIO) install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H" diff --git a/doc/header.html b/doc/header.html index 776a3e13..8642c69a 100644 --- a/doc/header.html +++ b/doc/header.html @@ -4,7 +4,7 @@ - + $projectname: $title diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp index 43fb9500..6e2a8892 100644 --- a/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.hpp @@ -28,45 +28,68 @@ Licence: namespace pFlow { - +/** + * Second order Adams-Bashforth integration method for solving ODE + * + * This is a one-step integration method and does not have prediction step. + * + */ class AdamsBashforth2 : public integration { protected: + /// dy at t-dt realx3PointField_D& dy1_; + /// Range policy for integration kernel (alias) using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, Kokkos::IndexType >; + public: - // type info + /// Type info TypeInfo("AdamsBashforth2"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsBashforth2( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + uniquePtr clone()const override + { + return makeUnique(*this); + } + + /// Destructor virtual ~AdamsBashforth2()=default; - // - add a virtual constructor + /// Add this to the virtual constructor table add_vCtor( integration, AdamsBashforth2, word); - //// - Methods - bool predict(real UNUSED(dt), realx3Vector_D& UNUSED(y), realx3Vector_D& UNUSED(dy)) override; + // - Methods - bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + bool predict( + real UNUSED(dt), + realx3Vector_D& UNUSED(y), + realx3Vector_D& UNUSED(dy)) override; + + bool correct( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; bool setInitialVals( const int32IndexContainer& newIndices, @@ -76,16 +99,21 @@ public: { return false; } + + /// Integrate on all points in the active range + bool intAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng); - uniquePtr clone()const override - { - return makeUnique(*this); - } - - 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 ); + bool intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ); }; diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp index e120941a..040a9e40 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.cpp @@ -31,26 +31,8 @@ pFlow::AdamsBashforth3::AdamsBashforth3 ) : integration(baseName, owner, pStruct, method), - /*dy1_( - owner.emplaceObject( - objectFile( - groupNames(baseName,"dy1"), - "", - objectFile::READ_IF_PRESENT, - objectFile::WRITE_ALWAYS), - pStruct, - zero3)), - dy2_( - owner.emplaceObject( - objectFile( - groupNames(baseName,"dy2"), - "", - objectFile::READ_IF_PRESENT, - objectFile::WRITE_ALWAYS), - pStruct, - zero3))*/ history_( - owner.emplaceObject( + owner.emplaceObject>( objectFile( groupNames(baseName,"AB3History"), "", diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp b/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp index 11c93d5c..6bd47b11 100644 --- a/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp +++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.hpp @@ -30,10 +30,10 @@ namespace pFlow struct AB3History { + TypeInfoNV("AB3History"); + realx3 dy1_={0,0,0}; realx3 dy2_={0,0,0}; - - TypeInfoNV("AB3History"); }; @@ -65,21 +65,21 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3) return str; } - +/** + * Third order Adams-Bashforth integration method for solving ODE + * + * This is a one-step integration method and does not have prediction step. + */ class AdamsBashforth3 : public integration { protected: + + /// Integration history + pointField& history_; - using HistoryFieldType = pointField; - //realx3PointField_D& dy1_; - - //realx3PointField_D& dy2_; - - // this is a device - HistoryFieldType& history_; - + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -91,23 +91,32 @@ public: // type info TypeInfo("AdamsBashforth3"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsBashforth3( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + + uniquePtr clone()const override + { + return makeUnique(*this); + } + /// Destructor virtual ~AdamsBashforth3()=default; - // - add a virtual constructor + /// Add this to the virtual constructor table add_vCtor( integration, AdamsBashforth3, word); - //// - Methods + // - Methods + bool predict( real UNUSED(dt), realx3Vector_D & UNUSED(y), @@ -126,18 +135,17 @@ public: return false; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool intAll(real dt, + /// 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, + bool intRange( + real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP ); diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp index 0a713ffb..547de93f 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp +++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.cpp @@ -32,7 +32,7 @@ pFlow::AdamsBashforth4::AdamsBashforth4 : integration(baseName, owner, pStruct, method), history_( - owner.emplaceObject( + owner.emplaceObject>( objectFile( groupNames(baseName,"AB4History"), "", diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp b/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp index 909d93ba..4f4a2914 100644 --- a/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp +++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.hpp @@ -30,11 +30,12 @@ namespace pFlow struct AB4History { + TypeInfoNV("AB4History"); + realx3 dy1_={0,0,0}; realx3 dy2_={0,0,0}; realx3 dy3_={0,0,0}; - - TypeInfoNV("AB4History"); + }; @@ -68,23 +69,21 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4) return str; } - +/** + * Fourth order Adams-Bashforth integration method for solving ODE + * + * This is a one-step integration method and does not have prediction step. + */ class AdamsBashforth4 : public integration { protected: - using HistoryFieldType = pointField; - //realx3PointField_D& dy1_; - - //realx3PointField_D& dy2_; - - - - // this is a device - HistoryFieldType& history_; + /// Integration history + pointField& history_; + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -93,32 +92,42 @@ protected: public: - // type info + /// Type info TypeInfo("AdamsBashforth4"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsBashforth4( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + + uniquePtr clone()const override + { + return makeUnique(*this); + } + /// Destructor virtual ~AdamsBashforth4()=default; - // - add a virtual constructor + /// Add a this to the virtual constructor table add_vCtor( integration, AdamsBashforth4, word); - //// - Methods + // - Methods + bool predict( real UNUSED(dt), realx3Vector_D & UNUSED(y), realx3Vector_D& UNUSED(dy)) override; - bool correct(real dt, + bool correct( + real dt, realx3Vector_D & y, realx3Vector_D& dy) override; @@ -131,18 +140,17 @@ public: return false; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool intAll(real dt, + /// 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, + bool intRange( + real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP ); diff --git a/src/Integration/AdamsBashforth5/AdamsBashforth5.cpp b/src/Integration/AdamsBashforth5/AdamsBashforth5.cpp index 327b5f2b..580f0c06 100644 --- a/src/Integration/AdamsBashforth5/AdamsBashforth5.cpp +++ b/src/Integration/AdamsBashforth5/AdamsBashforth5.cpp @@ -32,7 +32,7 @@ pFlow::AdamsBashforth5::AdamsBashforth5 : integration(baseName, owner, pStruct, method), history_( - owner.emplaceObject( + owner.emplaceObject>( objectFile( groupNames(baseName,"AB5History"), "", diff --git a/src/Integration/AdamsBashforth5/AdamsBashforth5.hpp b/src/Integration/AdamsBashforth5/AdamsBashforth5.hpp index ce8a5d19..99261588 100644 --- a/src/Integration/AdamsBashforth5/AdamsBashforth5.hpp +++ b/src/Integration/AdamsBashforth5/AdamsBashforth5.hpp @@ -30,12 +30,12 @@ namespace pFlow struct AB5History { + TypeInfoNV("AB5History"); + realx3 dy1_={0,0,0}; realx3 dy2_={0,0,0}; realx3 dy3_={0,0,0}; realx3 dy4_={0,0,0}; - - TypeInfoNV("AB5History"); }; @@ -71,18 +71,21 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5) return str; } - +/** + * Fifth order Adams-Bashforth integration method for solving ODE + * + * This is a one-step integration method and does not have prediction step. + */ class AdamsBashforth5 : public integration { protected: - using HistoryFieldType = pointField; - - // this is a device - HistoryFieldType& history_; + /// Integration history + pointField& history_; + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -91,32 +94,40 @@ protected: public: - // type info + /// Type info TypeInfo("AdamsBashforth5"); - //// - Constructors + // - Constructors + AdamsBashforth5( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + uniquePtr clone()const override + { + return makeUnique(*this); + } + virtual ~AdamsBashforth5()=default; - // - add a virtual constructor + /// Add this to the virtual constructor table add_vCtor( integration, AdamsBashforth5, word); - //// - Methods + // - Methods + bool predict( real UNUSED(dt), realx3Vector_D & UNUSED(y), realx3Vector_D& UNUSED(dy)) override; - bool correct(real dt, + bool correct( + real dt, realx3Vector_D & y, realx3Vector_D& dy) override; @@ -129,18 +140,17 @@ public: return false; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool intAll(real dt, + /// 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, + bool intRange( + real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP ); @@ -184,4 +194,4 @@ bool pFlow::AdamsBashforth5::intRange( } // pFlow -#endif //__integration_hpp__ +#endif // diff --git a/src/Integration/AdamsMoulton3/AdamsMoulton3.hpp b/src/Integration/AdamsMoulton3/AdamsMoulton3.hpp index 2e33ae11..349d62b2 100644 --- a/src/Integration/AdamsMoulton3/AdamsMoulton3.hpp +++ b/src/Integration/AdamsMoulton3/AdamsMoulton3.hpp @@ -28,19 +28,27 @@ Licence: namespace pFlow { - +/** + * Third order Adams-Moulton integration method for solving ODE + * + * This is a predictor-corrector integration method. + */ class AdamsMoulton3 : public integration { protected: + /// y at time t realx3PointField_D& y0_; + /// dy at time t realx3PointField_D& dy0_; + /// dy at time t-dt realx3PointField_D& dy1_; + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -48,29 +56,44 @@ protected: >; public: - // type info + /// Type info TypeInfo("AdamsMoulton3"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsMoulton3( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + uniquePtr clone()const override + { + return makeUnique(*this); + } + + /// Destructor virtual ~AdamsMoulton3()=default; - // - add a virtual constructor + /// Add this to the virtual constructor table add_vCtor( integration, AdamsMoulton3, word); - //// - Methods - bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + // - Methods - bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + bool predict( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; + + bool correct( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; bool setInitialVals( const int32IndexContainer& newIndices, @@ -81,20 +104,35 @@ public: return true; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// Prediction step on all points in the active range + bool predictAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng); + /// Prediction step on active points in the active range template - bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP); + bool predictRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP); - bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// 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 ); + bool intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP); }; @@ -104,7 +142,7 @@ bool AdamsMoulton3::predictRange( real dt, realx3Vector_D& y, realx3Vector_D& dy, - activeFunctor activeP ) + activeFunctor activeP) { auto d_dy = dy.deviceVectorAll(); auto d_y = y.deviceVectorAll(); @@ -141,7 +179,7 @@ bool pFlow::AdamsMoulton3::intRange( real dt, realx3Vector_D& y, realx3Vector_D& dy, - activeFunctor activeP ) + activeFunctor activeP) { auto d_dy = dy.deviceVectorAll(); @@ -176,4 +214,4 @@ bool pFlow::AdamsMoulton3::intRange( } // pFlow -#endif //__integration_hpp__ +#endif // diff --git a/src/Integration/AdamsMoulton4/AdamsMoulton4.hpp b/src/Integration/AdamsMoulton4/AdamsMoulton4.hpp index 1e3863a7..c0da0f62 100644 --- a/src/Integration/AdamsMoulton4/AdamsMoulton4.hpp +++ b/src/Integration/AdamsMoulton4/AdamsMoulton4.hpp @@ -28,21 +28,30 @@ Licence: namespace pFlow { - +/** + * Fourth order Adams-Moulton integration method for solving ODE + * + * This is a predictor-corrector integration method. + */ class AdamsMoulton4 : public integration { protected: + /// y at time t realx3PointField_D& y0_; + /// dy at time t realx3PointField_D& dy0_; + /// dy at time t-dt realx3PointField_D& dy1_; + /// dy at time t-2*dt realx3PointField_D& dy2_; + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -50,29 +59,44 @@ protected: >; public: - // type info + /// Type info TypeInfo("AdamsMoulton4"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsMoulton4( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + uniquePtr clone()const override + { + return makeUnique(*this); + } + + /// Destructor virtual ~AdamsMoulton4()=default; - // - add a virtual constructor + // Add this to the virtual constructor table add_vCtor( integration, AdamsMoulton4, word); - //// - Methods - bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + // - Methods - bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + bool predict( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; + + bool correct( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; bool setInitialVals( const int32IndexContainer& newIndices, @@ -83,20 +107,35 @@ public: return true; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// Prediction step on all points in the active range + bool predictAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng); + /// Prediction step on active points in the active range template - bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP); + bool predictRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP); - bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// 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 ); + bool intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ); }; @@ -182,4 +221,4 @@ bool pFlow::AdamsMoulton4::intRange( } // pFlow -#endif //__integration_hpp__ +#endif // diff --git a/src/Integration/AdamsMoulton5/AdamsMoulton5.hpp b/src/Integration/AdamsMoulton5/AdamsMoulton5.hpp index bcfe9035..2e381d4a 100644 --- a/src/Integration/AdamsMoulton5/AdamsMoulton5.hpp +++ b/src/Integration/AdamsMoulton5/AdamsMoulton5.hpp @@ -28,23 +28,33 @@ Licence: namespace pFlow { - +/** + * Fifth order Adams-Moulton integration method for solving ODE + * + * This is a predictor-corrector integration method. + */ class AdamsMoulton5 : public integration { protected: + /// y at time t realx3PointField_D& y0_; + /// dy at time t realx3PointField_D& dy0_; + /// dy at time t-dt realx3PointField_D& dy1_; + /// dy at time t-2*dt realx3PointField_D& dy2_; + /// dy at time t-3*dt realx3PointField_D& dy3_; + /// Range policy for integration kernel using rpIntegration = Kokkos::RangePolicy< DefaultExecutionSpace, Kokkos::Schedule, @@ -52,16 +62,24 @@ protected: >; public: - // type info + /// Type info TypeInfo("AdamsMoulton5"); - //// - Constructors + // - Constructors + + /// Construct from components AdamsMoulton5( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + + uniquePtr clone()const override + { + return makeUnique(*this); + } + /// Destructor virtual ~AdamsMoulton5()=default; // - add a virtual constructor @@ -71,10 +89,17 @@ public: word); - //// - Methods - bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + // - Methods + + bool predict( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; - bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + bool correct( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy) override; bool setInitialVals( const int32IndexContainer& newIndices, @@ -85,20 +110,35 @@ public: return true; } - uniquePtr clone()const override - { - return makeUnique(*this); - } - - bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// Prediction step on all points in the active range + bool predictAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng); + /// Prediction step on active points in the active range template - bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP); + bool predictRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP); - bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + /// 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 ); + bool intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ); }; diff --git a/src/Integration/integration/integration.hpp b/src/Integration/integration/integration.hpp index a076039b..3cca58c2 100644 --- a/src/Integration/integration/integration.hpp +++ b/src/Integration/integration/integration.hpp @@ -31,32 +31,71 @@ Licence: namespace pFlow { - +/** + * Base class for integrating the first order ODE (IVP) + * + * The ODE should be in the following form: + *\f[ + \frac{dy}{dt} = f(y,t) + \f] + * for example the equation of motion is in the following form: + *\f[ + m\frac{d\vec{v}}{dt} = m\vec{g} + \sum \vec{f_c}(\vec{v},t) + \f] + * + * The integration method can be either one-step or predictor-corrector type. + * + */ class integration { protected: - repository& owner_; + // - Protected data members - const word baseName_; + /// The owner repository that all fields are storred in + repository& owner_; - const pointStructure& pStruct_; + /// The base name for integration + const word baseName_; + + /// A reference to pointStructure + const pointStructure& pStruct_; public: - // type info + /// Type info TypeInfo("integration"); - //// - Constructors + + // - Constructors + + /// Construct from components integration( const word& baseName, repository& owner, const pointStructure& pStruct, const word& method); + /// Copy constructor + integration(const integration&) = default; + + /// Move constructor + integration(integration&&) = default; + + /// Copy assignment + integration& operator = (const integration&) = default; + + /// Move assignment + integration& operator = (integration&&) = default; + + /// Polymorphic copy/cloning + virtual + uniquePtr clone()const=0; + + /// Destructor virtual ~integration()=default; - // - add a virtual constructor + /// Add a virtual constructor create_vCtor( integration, word, @@ -67,35 +106,49 @@ public: (baseName, owner, pStruct, method) ); - //// - Methods + // - Methods + /// Const ref to pointStructure + inline const auto& pStruct()const { return pStruct_; } - virtual bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; - - virtual bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; - - virtual bool setInitialVals( - const int32IndexContainer& newIndices, - const realx3Vector& y) = 0; - - virtual bool needSetInitialVals()const = 0; - - virtual uniquePtr clone()const=0; - + /// Base name + inline const word& baseName()const { return baseName_; } + /// Ref to the owner repository + inline repository& owner() { return owner_; } + /// Prediction step in integration + virtual + bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; + + /// Correction/main integration step + virtual + bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0; + + /// Set the initial values for new indices + virtual + bool setInitialVals( + const int32IndexContainer& newIndices, + const realx3Vector& y) = 0; + + /// Check if the method requires any set initial vals + virtual + bool needSetInitialVals()const = 0; + + + /// Create the polymorphic object based on inputs static uniquePtr create( const word& baseName, @@ -103,7 +156,7 @@ public: const pointStructure& pStruct, const word& method); -}; +}; // integration } // pFlow diff --git a/src/Integration/integration/integrations.hpp b/src/Integration/integration/integrations.hpp deleted file mode 100644 index 428656c2..00000000 --- a/src/Integration/integration/integrations.hpp +++ /dev/null @@ -1,28 +0,0 @@ -/*------------------------------- 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 __integrations_hpp__ -#define __integrations_hpp__ - -#include "integration.hpp" -#include "AdamsBashforth2.hpp" -#include "AdamsBashforth3.hpp" - -#endif \ No newline at end of file diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp index 7e101b31..843a3c46 100644 --- a/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp +++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.hpp @@ -24,7 +24,7 @@ Licence: #include "Time.hpp" #include "pointFields.hpp" -#include "integrations.hpp" +#include "integration.hpp" #include "uniquePtr.hpp" namespace pFlow @@ -57,17 +57,6 @@ public: dynamicPointStructure(const dynamicPointStructure& ps) = default; - /*dynamicPointStructure(const dynamicPointStructure& ps): - pointStructure(ps), - time_(ps.time_), - integrationMethod_(ps.integrationMethod_), - velocity_(ps.velocity_), - integrationPos_(ps.integrationPos_->clone()), - integrationVel_(ps.integrationVel_->clone()) - { - - }*/ - // - no move construct dynamicPointStructure(dynamicPointStructure&&) = delete; diff --git a/src/demComponent/demComponent.hpp b/src/demComponent/demComponent.hpp index 14041057..d5ad1635 100644 --- a/src/demComponent/demComponent.hpp +++ b/src/demComponent/demComponent.hpp @@ -28,70 +28,132 @@ Licence: namespace pFlow { +/** + * A base class for every main component of DEM system. + * + * This class provides abstraction at a very high level for any solver/utility + * that forces the derived component (i.e. particles, geometry, and etc) to + * advance over time when iterate is called in time a loop. + * + */ class demComponent { protected: - word componentName_; + // - Protected data members + /// Name of the DEM component + word componentName_; - systemControl& control_; + /// Reference to systemControl + systemControl& control_; - - Timers timers_; + /// All timers (if any) of this component + Timers timers_; public: + /// Type info TypeInfo("demComponent"); - demComponent(const word& name, systemControl& control) - : - componentName_(name), - control_(control), - timers_(name, &control.timers()) - {} + // - Constructors + + /// construct from components + demComponent(const word& name, systemControl& control) + : + componentName_(name), + control_(control), + timers_(name, &control.timers()) + {} - virtual ~demComponent() = default; + /// No copy constructor + demComponent(const demComponent&) = delete; + /// No move constructor + demComponent(demComponent&&) = delete; - const auto& control()const - { - return control_; - } + /// No copy assignment + demComponent& operator = (const demComponent&) = delete; - auto& control() - { - return control_; - } + /// No move assignment + demComponent& operator =(demComponent&&) = delete; - inline - real dt()const - { - return control_.time().dt(); - } + /// destructor + virtual ~demComponent() = default; - inline - real currentTime()const - { - return control_.time().currentTime(); - } + + // - Member functions - auto& timers(){ - return timers_; - } + /// Const ref to systemControl + inline + const auto& control()const + { + return control_; + } - const auto& timers() const{ - return timers_; - } + /// Ref to systemControl + inline + auto& control() + { + return control_; + } + /// Time step of integration + inline + real dt()const + { + return control_.time().dt(); + } - virtual bool beforeIteration() = 0; + /// Current simulation time + inline + real currentTime()const + { + return control_.time().currentTime(); + } + /// Const ref to timers + inline + const auto& timers() const + { + return timers_; + } - virtual bool iterate() = 0; + /// Ref to timers + inline + auto& timers() + { + return timers_; + } + /// This is called before the start of time loop + virtual + bool beforeTimeLoop() + { + notImplementedFunction; + return false; + } - virtual bool afterIteration() = 0; + /// This is called in time loop, before iterate + virtual + bool beforeIteration() = 0; + + /// This is called in time loop. Perform the main calculations + /// when the component should evolve along time. + virtual + bool iterate() = 0; + + /// This is called in time loop, after iterate. + virtual + bool afterIteration() = 0; + + /// This is called after the time loop + virtual + bool afterTimeLoop() + { + notImplementedFunction; + return false; + } }; diff --git a/tutorials/sphereGranFlow/RotatingDrumWithBaffles/ReadMe.md b/tutorials/sphereGranFlow/RotatingDrumWithBaffles/ReadMe.md index d39e076e..8815cb83 100644 --- a/tutorials/sphereGranFlow/RotatingDrumWithBaffles/ReadMe.md +++ b/tutorials/sphereGranFlow/RotatingDrumWithBaffles/ReadMe.md @@ -1,5 +1,5 @@ # Problem Definition -The problem is to simulate a rotating drum with the diameter **0.24 m**, the length **0.1 m** and **6** Baffles, rotating at **15 rpm**. This drum is filled with **20000** Particles.The timestep for integration is **0.00001 s**. There are 2 types of Particles in this drum each are beining inserted during simulation to fill the drum. +The problem is to simulate a rotating drum with the diameter **0.24 m**, the length **0.1 m** and **6** Baffles, rotating at **15 rpm**. This drum is filled with **20000** Particles.The timestep for integration is **0.00001 s**. There are 2 types of Particles in this drum each are being inserted during simulation to fill the drum. * **12500** Particles with **4 mm** diameter, at the rate of 12500 particles/s for 1 sec. * **7500** Particles with **5mm** diameter, at the rate of 7500 particles/s for 1 sec. diff --git a/tutorials/sphereGranFlow/toteblender/ReadMe.md b/tutorials/sphereGranFlow/toteblender/ReadMe.md index 5a6d13b5..7b13e340 100644 --- a/tutorials/sphereGranFlow/toteblender/ReadMe.md +++ b/tutorials/sphereGranFlow/toteblender/ReadMe.md @@ -1,6 +1,6 @@ # Problem Definition -The problem is to simulate a double pedestal tote blender with the diameter **0.03 m** and **0.1 m** respectively, the length **0.3 m**, rotating at **28 rpm**. This blender is filled with **20000** Particles. The timestep for integration is **0.00001 s**. There is one type of Particle in this blender that are being inserted during simulation to fill the blender. -* **20000** particles with **4 mm** diameter, at the rate of 20000 particles/s for 1 sec. ŮŽAfter settling particles, this blender starts to rotate at t=**1s**. +The problem is to simulate a double pedestal tote blender (mixer) with the diameter **0.03 m** and **0.1 m** respectively, the length **0.3 m**, rotating at **28 rpm**. This blender is filled with **24000** particles. The timestep for integration is **0.00001 s**. There is one type of particle in this blender. Particles are positioned before start of simulation to fill the blender. +* **24000** particles with **5 mm** diameter are positioned, in order, and let to be settled under gravity. After settling particles, this blender starts to rotate at t=**1s**. @@ -8,238 +8,22 @@ The problem is to simulate a double pedestal tote blender with the diameter **0. a view of the tote-blender while rotating
- +
+
+ particles are colored according to their velocity +
# Setting up the Case -As it has been explained in the previous cases, the simulation case setup is based on text-based scripts. Here, the simulation case setup are sotred in two folders: `caseSetup`, `setting`. (see the above folders). Unlike the previous cases, this case does not have the `stl` file. and the geometry is described in the `geometryDict` file. +As it has been explained in the previous cases, the simulation case setup is based on text-based scripts. Here, the simulation case setup files are stored into two folders: `caseSetup`, `setting` (see the above folders). Unlike the previous cases, this case does not have the `stl` file and the surfaces are defined based on the built-in utilities in phasicFlow. See next the section for more information on how we can setup the geometry and its rotation. -## Defining particles -Then in the `caseSetup/sphereShape` the diameter and the material name of the particles are defined. -```C++ -// names of shapes -names (sphere1); -// diameter of shapes (m) -diameters (0.004); -// material names for shapes -materials (prop1); -``` -## Particle Insertion -In this case we have a region for ordering particles. These particles are placed in this blender. For example the script for the inserted particles is shown below. +## Geometry -
-in caseSetup/particleInsertion file -
+### Defining rotation axis +In file `settings/geometryDict` the information of rotating axis and speed of rotation are defined. The rotation of this blender starts at time=**0.5 s** and ends at time=**9.5 s**. -```C++ -// positions particles -positionParticles -{ -// ordered positioning - method positionOrdered; -// maximum number of particles in the simulation - maxNumberOfParticles 40000; -// perform initial sorting based on morton code? - mortonSorting Yes; -// cylinder for positioning particles - cylinder - { -// Coordinates of top cylinderRegion (m,m,m) - p1 (0.05 0.0 0.12); - p2 (0.05 0.0 0.22); -// radius of cylinder - radius 0.066; - } - - positionOrderedInfo - { -// minimum space between centers of particles - diameter 0.003; -// number of particles in the simulation - numPoints 20000; -// axis order for filling the space with particles - axisOrder (z y x); - } -} -``` - ## Interaction between particles - In `caseSetup/interaction` file, material names and properties and interaction parameters are defined: interaction between the particles of rotating drum. Since we are defining 1 material for simulation, the interaction matrix is 1x1 (interactions are symetric). -```C++ - // a list of materials names -materials (prop1); -// density of materials [kg/m3] -densities (1000.0); - -contactListType sortedContactList; - -model -{ - contactForceModel nonLinearNonLimited; - rollingFrictionModel normal; - /* - Property (prop1-prop1); - */ -// Young modulus [Pa] - Yeff (1.0e6); -// Shear modulus [Pa] - Geff (0.8e6); -// Poisson's ratio [-] - nu (0.25); -// coefficient of normal restitution - en (0.7); -// coefficient of tangential restitution - et (1.0); -// dynamic friction - mu (0.3); -// rolling friction - mur (0.1); - -} -``` -## Settings -### Geometry -In the `settings/geometryDict` file, the geometry and axis of rotation is defined for the drum. The geometry is composed of a cylinder inlet and outlet, cone shell top and down, a cylinder shell and enter and exit Gate. -```C++ -surfaces -{ - topGate - topGate - { - // type of wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.299); - // end point of cylinder axis - p2 (0.0 0.0 0.3); - // radius at p1 - radius1 0.03; - // radius at p2 - radius2 0.0001; - // material of wall - material solidProperty; - // motion component name - motion axisOfRotation; - } - - topCylinder - { - // type of the wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.28); - // end point of cylinder axis - p2 (0.0 0.0 0.3); - // radius at p1 - radius1 0.03; - // radius at p2 - radius2 0.03; - // number of divisions - resolution 36; - // material name of this wall - material prop1; - // motion component name - motion axisOfRotation; - } - - coneShelltop - { - // type of the wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.2); - // end point of cylinder axis - p2 (0.0 0.0 0.28); - // radius at p1 - radius1 0.1; - // radius at p2 - radius2 0.03; - // number of divisions - resolution 36; - // material name of this wall - material prop1; - // motion component name - motion axisOfRotation; - } - - cylinderShell - { - // type of the wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.1); - // end point of cylinder axis - p2 (0.0 0.0 0.2); - // radius at p1 - radius1 0.1; - // radius at p2 - radius2 0.1; - // number of divisions - resolution 36; - // material name of this wall - material prop1; - // motion component name - motion axisOfRotation; - } - - coneShelldown - { - // type of the wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.02); - // end point of cylinder axis - p2 (0.0 0.0 0.1); - // radius at p1 - radius1 0.03; - // radius at p2 - radius2 0.1; - // number of divisions - resolution 36; - // material name of this wall - material prop1; - // motion component name - motion axisOfRotation; - } - /* - This is a plane wall at the exit of silo - */ - - bottomCylinder - { - // type of the wall - type cylinderWall; - // begin point of cylinder axis - p1 (0.0 0.0 0.0); - // end point of cylinder axis - p2 (0.0 0.0 0.02); - // radius at p1 - radius1 0.03; - // radius at p2 - radius2 0.03; - // number of divisions - resolution 36; - // material name of this wall - material prop1; - // motion component name - motion axisOfRotation; - } - exitGate - { - type planeWall; - p1 (-0.05 -0.05 0); - p2 (-0.05 0.05 0); - p3 ( 0.05 0.05 0); - p4 (0.05 -0.05 0); - material prop1; - motion axisOfRotation; - } - -} -``` -### Rotating Axis Info -In this part of `geometryDict` the information of rotating axis and speed of rotation are defined. Unlike the previous cases, the rotation of this blender starts at time=**0 s**. ```C++ // information for rotatingAxisMotion motion model rotatingAxisMotionInfo @@ -247,19 +31,326 @@ rotatingAxisMotionInfo axisOfRotation { p1 (-0.1 0.0 0.15); // first point for the axis of rotation - p2 (0.1 0.0 0.15); // second point for the axis of rotation + p2 ( 0.1 0.0 0.15); // second point for the axis of rotation + omega 1.5708; // rotation speed ==> 15 rad/s - // Start time of Geometry Rotating (s) - startTime 1; - // End time of Geometry Rotating (s) + + // Start time of Geometry Rotating (s) + startTime 0.5; + + // End time of Geometry Rotating (s) endTime 9.5; } } ``` -## Performing Simulation + + +### Surfaces +In `settings/geometryDict` file, the surfaces and motion component of each surface are defined to form a rotating tote-blender. The geometry is composed of top and bottom cylinders, top and bottom cones, a cylindrical shell and top and bottom Gates. + +```C++ +surfaces +{ + + topGate + { + // type of wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.3); + + // end point of cylinder axis + p2 (0.0 0.0 0.301); + + // radius at p1 + radius1 0.03; + + // radius at p2 + radius2 0.0001; + + // material of wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + topCylinder + { + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.28); + + // end point of cylinder axis + p2 (0.0 0.0 0.3); + + // radius at p1 + radius1 0.03; + + // radius at p2 + radius2 0.03; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + coneShelltop + { + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.2); + + // end point of cylinder axis + p2 (0.0 0.0 0.28); + + // radius at p1 + radius1 0.1; + + // radius at p2 + radius2 0.03; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + cylinderShell + { + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.1); + + // end point of cylinder axis + p2 (0.0 0.0 0.2); + + // radius at p1 + radius1 0.1; + + // radius at p2 + radius2 0.1; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + coneShelldown + { + + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.02); + + // end point of cylinder axis + p2 (0.0 0.0 0.1); + + // radius at p1 + radius1 0.03; + + // radius at p2 + radius2 0.1; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + bottomCylinder + { + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 0.0); + + // end point of cylinder axis + p2 (0.0 0.0 0.02); + + // radius at p1 + radius1 0.03; + + // radius at p2 + radius2 0.03; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + + exitGate + { + + // type of the wall + type cylinderWall; + + // begin point of cylinder axis + p1 (0.0 0.0 -0.001); + + // end point of cylinder axis + p2 (0.0 0.0 0.0); + + // radius at p1 + radius1 0.03; + + // radius at p2 + radius2 0.0001; + + // number of divisions + resolution 36; + + // material name of this wall + material solidProperty; + + // motion component name + motion axisOfRotation; + } + +} +``` + +## Defining particles +### Diameter and material of spheres +In the `caseSetup/sphereShape` the diameter and the material name of the particles are defined. + +
+in caseSetup/sphereShape file +
+ +```C++ +// name of shapes +names (sphere1); + +// diameter of shapes (m) +diameters (0.005); + +// material name for shapes +materials (solidProperty); +``` +### Particle positioning before start of simulation +Particles are positioned before the start of simulation. The positioning can be ordered or random. Here we use ordered positioning. 24000 particles are positioned in a cylinderical region inside the tote-blender. + +
+in settings/particlesDict file +
+ +```C++ +// positions particles +positionParticles +{ + // ordered positioning + method positionOrdered; + + // maximum number of particles in the simulation + maxNumberOfParticles 25001; + + // perform initial sorting based on morton code? + mortonSorting Yes; + + // cylinderical region for positioning particles + cylinder + { + p1 (0.0 0.0 0.09); + p2 (0.0 0.0 0.21); + radius 0.09; + } + + positionOrderedInfo + { + // minimum space between centers of particles + diameter 0.005; + + // number of particles in the simulation + numPoints 24000; + + // axis order for filling the space with particles + axisOrder (x y z); + } +} +``` + + ## Interaction between particles + In `caseSetup/interaction` file, material names and properties and interaction parameters are defined. Since we are defining 1 material type in the simulation, the interaction matrix is 1x1 (interactions are symmetric). +```C++ + // a list of materials names +materials (solidProperty); + +// density of materials [kg/m3] +densities (1000.0); + +contactListType sortedContactList; + +model +{ + contactForceModel nonLinearNonLimited; + + rollingFrictionModel normal; + + /* + Property (solidProperty-solidProperty); + */ + + // Young modulus [Pa] + Yeff (1.0e6); + + // Shear modulus [Pa] + Geff (0.8e6); + + // Poisson's ratio [-] + nu (0.25); + + // coefficient of normal restitution + en (0.7); + + // coefficient of tangential restitution + et (1.0); + + // dynamic friction + mu (0.3); + + // rolling friction + mur (0.1); +} +``` + +# Performing Simulation and previewing the results To perform simulations, enter the following commands one after another in the terminal. Enter `$ particlesPhasicFlow` command to create the initial fields for particles. -Enter `$ geometryPhasicFlow` command to create the Geometry. +Enter `$ geometryPhasicFlow` command to create the geometry. At last, enter `$ sphereGranFlow` command to start the simulation. -After finishing the simulation, you can use `$ pFlowtoVTK` to convert the results into vtk format storred in ./VTK folder. \ No newline at end of file +After finishing the simulation, you can use `$ pFlowtoVTK` to convert the results into vtk format stored in ./VTK folder.