Phi (volume fraction) returned from getDepthProfiles were far too low

Hey all,

At this point, I am perplexed to see volume fraction (phi) from getDepthProfile returned from each cell (10 cells total) so incredibly low (roughly 2%) when the simulation of the cylindrical model shows extensive compaction of the spheres and should be much higher. Here is a few snippets of code to show where the work process has gone. I believe that this has something to do with the direction or perhaps the void volume fraction? Maybe the ram motion (for compaction) triggers this anomaly?

Start of code:

from yade import ymport, plot, utils
import time
import numpy as np

#Constants

utils.readParamsFromTable(descriptionIn = 'ram',
  frIn = 0.5, enIn=0.01, etIn=0.01, tcIn=0.0001, # Frictionangle (rad), strength-parameter
  rhoIn = 2300.0, # Density
  dumpVTKIn = 1000, # Periods of dumps
  ram_vel = 1E-1, ram_mov = 0.1E-3 # Ram Velocity and Distance
)

#Sweep Parameters using Table
from yade.params.table import *

import shutil

try:
      shutil.rmtree('step2')
except OSError:
      pass
os.mkdir('step2')

#Create Folder

folderNameBase = 'step2/' + str(descriptionIn)
folderName = folderNameBase
os.mkdir(folderNameBase)

#Initialize
o = Omega()
o.dt = 0.05*1E-7

#Read in previous step
o.load('save1.xml')

#remember start time
starttime = o.time #start time
stoptime = starttime+ram_mov/ram_vel #stop time
iterstop=o.iter+int((ram_mov/ram_vel)/o.dt)
print('start,stop,steps=',starttime,stoptime,iterstop)

#Define Next Scene
#o.addScene()

#Build Geometry
ram=o.bodies.append(ymport.stl('coax_ram.STL'))

o.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],sortThenCollide=False,verletDist=-0.5,label='collider'),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(damping=0.5,gravity=Vector3(0,-90.81,0), label='newt'),

   #Record output to VTK Files, Spheres,Interactions,Facets

  VTKRecorder(iterPeriod=dumpVTKIn,fileName=folderName+'/spheres-',
        recorders=['spheres','velocity','colors','intr','ids','mask','materialId','stress'], label='VTK_sphere'),
    VTKRecorder(iterPeriod=dumpVTKIn,fileName=folderName+'/facetB-',
        recorders=['facets','colors'],label='VTK_OD'),

#FIXME Snapshot Taker

# SnapshotEngine(iterPeriod=dumpVTKIn,fileBase=folderName+'/spheres-',viewNo=0,label='snapshooter'),

   #Execute Python Commands
    PyRunner(iterPeriod=10000,initRun=True,command='stop_condition()'),
        PyRunner(iterPeriod=iterstop,command='o.save(\'save2.xml\')'),
    PyRunner(iterPeriod=10000,initRun=True,command='addPlotPhi()')

#Enable Track Energy - Can be used to determine if sim is done
o.trackEnergy=True

#Run Yade Simulation
#o.stopAtIter=2150001 #Maximum Iterations
o.stopAtIter=iterstop+1
#o.stopAtTime=stoptime+1E-8
o.run()

def stop_condition():
    o.dt=0.5*utils.PWaveTimeStep() #auto timesetup
    phiavg=calcPhi()
    print('iter,time,dt,iter/sec,phi=',o.iter,o.iter*o.dt,o.dt,o.speed,phiavg) #print to screen iteration number

def calcPhi():
    CN=utils.getDepthProfiles(6.28E-8,10,2.0E-4,0.00,False,1E-6,2)
    return np.mean(CN[0])

def addPlotData():
    plot.addData(i=o.iter,unbalanced=unbalancedForce(),**o.energy)

--------------------------------------------------------------------------------------------------------------------------------
######Now this is the part of the code that I think may be why the phiPart returned is so low:
----------------------------------------------------------------------------------------------------------------------------------

def addPlotPhi():
    #getDepthProfiles(double volume, int nCell, double dz, double zRef, bool activateCond, double radiusPy, int direction)
    CN=utils.getDepthProfiles(6.28E-8,10,2.0E-4,0.00,False,1E-6,2)
    phiavg=np.mean(CN[0])
    zarray=np.linspace(0,2e-3,10,endpoint=True)
  fname='ramvolfrac_'+str(O.iter)+'.csv'
    phidata = np.column_stack((zarray, CN[0]))
    np.savetxt(fname, phidata, delimiter=',')

------------------------------------------------------------------------------------------------------------------------------------

As for the parameters, they are listed in meters for inner and outer cylinder (outer_cylinder_radius=0.007m, inner_cylinder_radius=0.003m) with a height of 0.0015 m. I would like to have 10 cells (as seen in the addPlotPhi) and hopefully a realistic approximation of solid volume fraction. Thank you for any help!

-Kevin