global damping is added to the code

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<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+			d_y[i] = damping*(d_dy[i] + dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(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<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+				d_y[i] =  damping*(d_y[i] + dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(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);

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, 

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<real>(23.0 / 12.0) * d_dy[i] 
+			d_y[i] = damping*( d_y[i]+ dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
 						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
-						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+						+ static_cast<real>(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<real>(23.0 / 12.0) * d_dy[i] 
+				d_y[i] = damping * (d_y[i] + dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
 						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
-						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+						+ static_cast<real>(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);

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, 

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<real>(55.0 / 24.0) * d_dy[i]
 				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
 				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
 				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
-				);
+				));
 			d_dy3[i] = d_dy2[i];
 			d_dy2[i] = d_dy1[i];
 			d_dy1[i] = d_dy[i];
@@ -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<real>(55.0 / 24.0) * d_dy[i]
 				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
 				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
 				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
-				);
+				));
 				d_dy3[i] = d_dy2[i];
 				d_dy2[i] = d_dy1[i];
 				d_dy1[i] = d_dy[i];
@@ -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);

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, 

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;

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

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 ..."<<END_REPORT;
+	velDamping_ = makeUnique<globalDamping>(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;	
 }

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<integration>  integrationVel_ = nullptr;
 
+	uniquePtr<globalDamping> velDamping_ = nullptr;
+
 	Timer 					velocityUpdateTimer_;
 
 	/// @brief integration method for velocity and position
@@ -89,6 +93,8 @@ public:
 	/// 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);
 

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<real>("dampingFactor", static_cast<real>(1.0));
+
+    dampingFactor_ = max( dampingFactor_ , static_cast<real>(0.01));
+
+    performDamping_ = !equal(dampingFactor_, static_cast<real>(1.0));
+
+    if( performDamping_ )
+		REPORT(2)<<"Global damping "<<Yellow_Text("is active")<<
+		" and damping factor is "<<dampingFactor_<<END_REPORT;
+	else
+		REPORT(2)<<"Global damping "<<Yellow_Text("is not active")<<"."<<END_REPORT;
+
+}
+
+
+/*void pFlow::globalDamping::applyDamping
+(
+    const timeInfo& ti, 
+    realx3PointField_D& velocity
+)
+{
+	if(!performDamping_) return;
+	if(!timeControl_.timeEvent(ti.iter(), ti.t(), ti.dt()) )return;
+
+	auto d_v = velocity.deviceView();
+	auto activeRng = velocity.activeRange();
+	auto dmpng = dampingFactor_;
+
+	Kokkos::parallel_for(
+		"globalDamping::applyDamping",
+		deviceRPolicyStatic(activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			d_v[i] *= dmpng;
+		});
+	Kokkos::fence();
+	//REPORT(1)<<"Applied global damping "<<END_REPORT;
+}*/
+
+pFlow::real pFlow::globalDamping::dampingFactor(const timeInfo& ti)const
+{
+	if(!performDamping_) return 1.0;
+	if(!timeControl_.timeEvent(ti.iter(), ti.t(), ti.dt()) )return 1.0;
+	return dampingFactor_;
+}

src/Particles/globalDamping/globalDamping.hpp

@@ -0,0 +1,65 @@
+/*------------------------------- 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 __globalDamping_hpp__
+#define __globalDamping_hpp__
+
+#include "systemControl.hpp"
+#include "pointFields.hpp"
+#include "baseTimeControl.hpp"
+
+namespace pFlow
+{
+
+
+class globalDamping
+{
+private:
+
+    bool            performDamping_ = false;
+
+    real            dampingFactor_;
+
+    baseTimeControl timeControl_;
+
+
+public:
+
+
+    globalDamping(const systemControl& control);
+
+    ~globalDamping()=default;
+
+    //void applyDamping( const timeInfo& ti, realx3PointField_D& velocity);
+
+    bool performDamping()const
+    {
+        return performDamping_;
+    }
+
+    real dampingFactor(const timeInfo& ti)const;
+    
+};
+
+
+}
+
+
+#endif 

src/phasicFlow/repository/Time/baseTimeControl.cpp

@@ -65,6 +65,50 @@ pFlow::baseTimeControl::baseTimeControl
     
 }
 
+pFlow::baseTimeControl::baseTimeControl
+(
+	const dictionary& dict,
+	const real 	defInterval,
+	const word& intervalPrefix,
+	const real 	defStartTime
+)
+:
+	intervalPrefix_(intervalPrefix)
+{
+	auto tControl = dict.getVal<word>("timeControl");
+	if(tControl == "timeStep")
+		isTimeStep_ = true;
+	else if( tControl == "runTime" || 
+			 tControl == "simulationTime")
+		isTimeStep_ = false;
+	else
+	{
+		fatalErrorInFunction<<
+		"bad input for intervalControl in dictionary "<< dict.globalName()<<endl<<
+		"valid inputs are "<< 
+		wordList({"timeStep", "runTime", "simulationTime"})<<endl;
+		fatalExit;
+	}
+
+	word intervalWord = intervalPrefix.size()==0? word("interval"): intervalPrefix+"Interval";
+
+	if(!isTimeStep_)
+	{
+		auto startTime = (dict.getValOrSet<real>("startTime", defStartTime));
+    	auto endTime = (dict.getValOrSet<real>("endTime", largeValue));
+		auto interval = dict.getValOrSet<real>(intervalWord, defInterval);
+		rRange_ = realStridedRange(startTime, endTime, interval);
+
+	}
+	else
+	{
+		auto startTime = (dict.getValOrSet<int32>("startTime", 0));
+    	auto endTime = (dict.getValOrSet<int32>("endTime", largestPosInt32));
+		auto interval = dict.getValOrSet<int32>(intervalWord, 1);
+		iRange_ = int32StridedRagne(startTime, endTime, interval);
+	}
+}
+
 pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix)
 :
 	isTimeStep_(true),

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,