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,