How to output a list throughout timestep

Hi all,

I am simulating a direct shear experiment with a granular layer.
I would like to see how the force in the X direction (O.force(b).[1]) transmit within the granular layer during initial shearing.
I define a weighted average function called "getForce" to calculate the weighted average force of each particle at a certain Y position (Ystep).
Then, I define another function called "forceExport" to save the weighted force and the corresponded depth as a list with "np.savetxt".
I would like to output the list every 1000 iteration so I use Pyrunner to do this.
However, I cannot think of a way to save the list with a filename that is corresponded to the iteration period because it would overwrite the previous file.
Is there another way to save a list throughout the iteration period?
Thanks!

Here's part of my code:
###
O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(krot=roll_stiff,eta=roll_fric)],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=TSSC,defaultDt=-1)
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ,PyRunner(command='dataRecorder()',iterPeriod=1000,label='recData',dead=True)
 ,PyRunner(command="forceExport()",iterPeriod=1000,label='forceExp',dead=True)
 ]

def weight_function_force(y,y1,coeff):
    return(math.exp(-(((y-y1)**2)/(2*(coeff**2)))))

def getForce(y1):
  weight=0
  weightedForce=0
  totalWeightedForce=0
  totalWeight=0
  for i in range(2,len(O.bodies)):
    weight = weight_function_force(O.bodies[i].state.pos[1], y1, coeff)
    weightedForce = O.forces.f(i).[0] * weight
    totalWeightedForce = totalWeightedForce + weightedForce
    totalWeight = totalWeight + weight
  return (totalWeightedForce/totalWeight)

def forceExport():
    dh = (O.bodies[0].state.pos[1] - O.bodies[1].state.pos[1]) / Ystep
    list_force = []
    for i in range(Ystep):
        y1 = dh*i + O.bodies[1].state.pos[1]
        list_force.append([getForce(y1),y1])
    np.savetxt(OUT+'-Force',list_force, delimiter=',')
###

Cheers,
Chien-Cheng