yade-users team mailing list archive

[Raphaël Maurin <question253112@xxxxxxxxxxxxxxxxxxxxx>]

New question #253112 on Yade:
https://answers.launchpad.net/yade/+question/253112

Hi all, 

I modified the source code of yade to have an engine applying drag and lift force on each particle at each time step in function of a given average velocity field. 
The fluid flow is turbulent so that I determine every dt_f also randomly a fluid velocity fluctuation associated to each particle. (dt_f~100 dt)
To evaluate the drag and lift, I need then to pass the value of the fluctuation associated to each particles from the python script to the C++ engine. This takes the form of a vector of size len(O.bodies) and should then be given to the C++ engines every dt_f

However, when the number of particle is high (~50 000), the simulation stops after about 60s and it is written :
"processus arrêté" (processus stopped in french). 
This is most probably due to the fact the process is taking too much memory. For example a simulation of 83000 particles after 60s of virtual time simulated takes 14GB of memory and 10GB of SWAP. 

I made a test script and I found that the increase in memory consumption with time is due to the number of time I am passing the vector with the fluctuation associated to the particle from the python script to the C++ engine. 

I have possibilities to pass a smaller vector to the C++ engine, or to evaluate the turbulent fluctuation inside the C++ engine, however it does not seem normal to me to have memory leakage when passing a variable from python to C++.  

Any idea of what it is due to ? Is there a way to fix this problem ? 

Thanks for your help !

Raphaël

Yade version : 2014-06-29.git-de4c01a
linux version : Ubuntu 12.04
Hereafter the modification of the C++ code I made. 

ForceEngine.cpp : 

void HydroForceEngine::action(){
	FOREACH(Body::id_t id, ids){
		Body* b=Body::byId(id,scene).get();
		if (!b) continue;
		if (!(scene->bodies->exists(id))) continue;
		const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
		if (sphere){
			Vector3r posSphere = b->state->pos;//position vector of the sphere
			int p = floor((posSphere[2]-zRef)/deltaZ); //cell number in which the particle is
			if ((p<nCell)&&(p>0)) {
				Vector3r liftForce = Vector3r::Zero();
				Vector3r dragForce = Vector3r::Zero();

				Vector3r vFluid(vxFluid[p]+vxFluct[id],0.0,vzFluct[id]); //fluid velocity at this point (including fluctuations)
				Vector3r vPart = b->state->vel;
				Vector3r vRel = vFluid - vPart;
				//Drag force calculation
				Real Rep = vRel.norm()*sphere->radius*2*rhoFluid/viscoDyn; 
				Real A = sphere->radius*sphere->radius*Mathr::PI;	//Crossection of the sphere
				if (vRel.norm()!=0.0) {
					Real hindranceF = pow(1-phiPart[p],-expoRZ); //hindrance function
					Real Cd = (0.44 + 24.4/Rep)*hindranceF; //drag coefficient
					dragForce = 0.5*rhoFluid*A*Cd*vRel.squaredNorm()*vRel.normalized();
				}
				//lift force calculation due to difference of pressure (Saffman lift)
				int intRadius = floor(sphere->radius/deltaZ);
				if ((p+intRadius<nCell)&&(p-intRadius>0)&&(lift==true)) {
					Real vRelTop = vxFluid[p+intRadius] - vPart[0]; // relative velocity of the fluid wrt the particle at the top of the particle
					Real vRelBottom = vxFluid[p-intRadius] - vPart[0]; // same at the bottom
					liftForce[2] = 0.5*rhoFluid*A*Cl*(vRelTop*vRelTop-vRelBottom*vRelBottom);
				}
				//Archimedes force calculation
				Vector3r archimedesForce = -4.0/3.0*Mathr::PI*sphere->radius*sphere->radius*sphere->radius*rhoFluid*gravity;
				//add the force to the particle
				scene->forces.addForce(id,dragForce+liftForce+archimedesForce);		
			}
		}
	}
}



ForceEngine.hpp : 

class HydroForceEngine: public PartialEngine{
	public:
		virtual void action();
	YADE_CLASS_BASE_DOC_ATTRS(HydroForceEngine,PartialEngine,"Apply drag and lift force (and Archimedes force) due to a fluid flow vector (1D) to each sphere. The applied force reads\n\n.. math:: F_{d}=\\frac{1}{2} C_d A\\rho|\\vec{v_f - v}| vec{v_f - v} \n\n where $\\rho$ is the medium density (:yref:`density<HydroForceEngine.rhoFluid>`), $v$ is particle's velocity,  $v_f$ is the velocity of the fluid at the particle center,  $A$ is particle projected area (disc), $C_d$ is the drag coefficient. The formulation of the drag coefficient depends on the local particle reynolds number and the solid volume fraction. The formulation of the drag is Dallavalle with a correction of Richardson-Zaki to take into account the hindrance effect. This law is classical in sediment transport.\n The formulation of the lift is taken from Wiberg and Smith 1985 and is such as : \n\n.. math:: F_{L}=\\frac{1}{2} C_L A\\rho((v_f - v)^2{top} - (v_f - v)^2{bottom}) \n\n Where the subscript top and bottom means evaluated at the top (respectively the bottom) of the sphere considered. This formulation of the lift is a reformulation of the so-called saffman lift which is due to the difference of pressure inside a shear flow.",
		((Real,rhoFluid,1000,,"Density of the medium (fluid or air), by default - density of water"))
		((Real,viscoDyn,1e-3,,"Dynamic viscosity of the fluid"))
		((Real,zRef,,,"Position of the reference point which correspond to the first value of the fluid velocity"))
		((Real,nCell,,,"Size of the vector of the fluid velocity"))
		((Real,deltaZ,,,"Size of the discretization/the cell along z"))
		((Real,expoRZ,3.1,,"Value of the Richardson-Zaki exponent, for the correction due to hindrance"))
                ((Real,lift,true,,"Option to activate or not the evaluation of the lift"))
		((Real,Cl,0.2,,"Value of the lift coefficient"))
                ((Vector3r,gravity,,,"Gravity vector"))
		((vector<Real>,vxFluid,,,"Streamwise fluid velocity profile in function of the altitude"))
		((vector<Real>,phiPart,,,"Solid volume fraction profile in function of the altitude"))
		((vector<Real>,vxFluct,,,"Vector containing the value of the fluctuation of the streamwise fluid velocity at the position of the particles considered. The ith component of this vector correspond to the fluctuation of fluid velocity along x at the position of particle of id i."))
		((vector<Real>,vzFluct,,,"Same as :yref:`vxFluct<DragEngineNEW.vxFluct>` but for the vertical direction z."))
	);
};
REGISTER_SERIALIZABLE(HydroForceEngine);


-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.