memory leakage when passing a variable from python to C++

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);