Accumulation effect of heat over time steps

Hi Jan,

Thanks for your code and prompt reply! It works well.

So if I want to export these accumulate values and visualize in the Paraview, what could I do?
I would be grateful if you could provide me more instructions.

Cheers,
Chien-Cheng

Updated MWE:
### assuming Python3, Python2 can be adjusted
from __future__ import print_function
from yade import export

class SomeCalculator:
   def __init__(self):
      self.cumulativeData = dict((a,0) for a in range(0,3))
      self.currentData = dict(self.cumulativeData)
   def collect(self):
      """Collect data from simulation"""
      # acumulate cumulative data
      for k,v in self.currentData.items():
         self.cumulativeData[k] += v
      # compute new current data
      self.currentData = dict((a,self.getDataFromOneParticle(a)) for a in range(0,3))
   def getDataFromOneParticle(self,a):
      b = O.bodies[a]
      ret = 0
      for i in b.intrs():
         if not i.isReal: continue
         ret += i.geom.penetrationDepth
      return ret
   def pprint(self):
      print()
      print("current step")
      print("============")
      for k,v in self.currentData.items():
         print(k,v)
      print("cumulative")
      print("==========")
      for k,v in self.cumulativeData.items():
         print(k,v)
   def getCurrent(self,b):
      return self.currentData[b.id]
   def getCumulative(self,b):
      return self.cumulativeData[b.id]

O.engines = [
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()],
   ),
   NewtonIntegrator(),
   PyRunner(iterPeriod=1,command="someCalculator.collect()"),
   PyRunner(iterPeriod=1,command="someCalculator.pprint()"),
   PyRunner(iterPeriod=1,command="vtkExport()"),
]

s1,s2,s3 = (
   sphere((0,0,0),1),
   sphere((1,0,0),1),
   sphere((1,1,0),1),
)
O.bodies.append((s1,s2,s3))

vtk = export.VTKExporter("someData")
def vtkExport():
   vtk.exportSpheres(what=dict(current="someCalculator.getCurrent(b)",cumulative="someCalculator.getCumulative(b)"))

# currently needs to be after bodies are created, could be adjusted..
someCalculator = SomeCalculator()
import builtins
builtins.someCalculator = someCalculator # such that VTKExporter has acces to it

O.dt = 0.003
O.run(10,True)
###

cheers
Jan

Hi Jan,

Many thanks again!

This time I got some errors as below. The cumulative data can still be printed out but the exported vtk file cannot be opened in the Paraview.

---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
/home/hung/YadeNov2019/bins/yade-2019-11-18.git-b1a2767 in <module>()

/home/hung/YadeNov2019/bins/yade-2019-11-18.git-b1a2767 in vtkExport()
    223 vtk = export.VTKExporter("someData")
    224 def vtkExport():
--> 225 vtk.exportSpheres(what=dict(current="someCalculator.getCurrent(b)",cumulative="someCalculator.getCumulative(b)"))
    226
    227 import builtins

/home/hung/YadeNov2019/build/lib/x86_64-linux-gnu/yade-2019-11-18.git-b1a2767/py/yade/export.py in exportSpheres(self, ids, what, comment, numLabel, useRef)
    427 # write additional data from 'what' param
    428 for name,command in what.items(): # for each name...
--> 429 test = eval(command) # ... eval one example to see what type (float, Vector3, Matrix3) the result is ...
    430 # ... and write appropriate header line and loop over all bodies and write appropriate vtk line(s)
    431 if isinstance(test,Matrix3):

/home/hung/YadeNov2019/build/lib/x86_64-linux-gnu/yade-2019-11-18.git-b1a2767/py/yade/export.py in <module>()

/home/hung/YadeNov2019/bins/yade-2019-11-18.git-b1a2767 in getCumulative(self, b)
    127 return self.currentData[b.id]
    128 def getCumulative(self,b):
--> 129 return self.cumulativeData[b.id]
    130
    131

KeyError: 1666
---------------------------------------------------------------------------

Do you have any ideas about this?

Cheers,
Chien-Cheng

> the exported vtk file cannot be opened in the Paraview.

because of the errors, no output is generated

> Do you have any ideas about this?

you try to access data, which are not present, a complete script would be needed to answer why..

cheers
Jan

Dear Jan,

I simplified my script and hope that you could have a look.

#######
from __future__ import print_function
import builtins
from yade import pack,plot,export
import numpy as np
import math

sp1=pack.SpherePack()
sp2=pack.SpherePack()
sp3=pack.SpherePack()

O.periodic=True

# dimensions of sample
RADIUS=0.025
a=15
b=1
c=1
length=a*(2*RADIUS)
height=length/b
width=length/c
thickness=RADIUS
psdSizes=[.0456,.05,.0544]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=0
spFRIC=26.6

# boundary conditions
PI=1.e5
SN=5.e6
RATE_shear=1

# simulation control
DAMPSHEAR=0.

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=radians(wallFRIC),label='boxMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=radians(spFRIC),label='sphereMat'))

upBox = utils.box(center=(length/2,2*height+thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

sp1.makeCloud((0,height+3*RADIUS,width),(length,2*height-3*RADIUS,2*width), psdSizes =psdSizes, psdCumm =psdCumm, periodic=True, seed =1)
sp2.makeCloud((0,height+RADIUS,width),(length,height+RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1)
sp3.makeCloud((0,2*height-RADIUS,width),(length,2*height-RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1)

sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp1])
bottomLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='sphereMat') for s in sp2])
topLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='sphereMat') for s in sp3])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

class SomeCalculator:
    def __init__(self):
        self.cumulativeData = dict((a,0) for a in range(2,10))
        self.currentData = dict(self.cumulativeData)
    def collect(self):
        """Collect data from simulation"""
        # acumulate cumulative data
        for key,valve in self.currentData.items():
           self.cumulativeData[key] += valve
        # compute new current data
        self.currentData = dict((a,self.getDataFromOneParticle(a)) for a in range(2,10))
    def getDataFromOneParticle(self,a):
        b = O.bodies[a]
        totalHeat = 0
        for i in b.intrs():
            O.interactions.erase(i.id1,0)
            O.interactions.erase(i.id1,1)
            O.interactions.erase(0,i.id2)
            O.interactions.erase(1,i.id2)
            if not i.isReal: continue
            penetrationDepth = i.geom.penetrationDepth
            radius1 = i.geom.refR1
            radius2 = i.geom.refR2
            effectiveRadius = (radius1*radius2)/(radius1+radius2)
            contactRadius = math.sqrt(penetrationDepth*effectiveRadius)
            contactArea = np.pi*(contactRadius**2)
            relativeVelocity = i.geom.shearInc.norm()/O.dt
            generatedHeat = (i.phys.shearForce.norm()/contactArea)/(i.phys.normalForce.norm()/contactArea) * i.phys.normalForce.norm() * relativeVelocity
            totalHeat = totalHeat + generatedHeat
        return totalHeat
    def pprint(self):
        print()
        print("current step")
        print("============")
        for key,valve in self.currentData.items():
            print(key,valve)
        print("cumulative")
        print("==========")
        for key,valve in self.cumulativeData.items():
            print(key,valve)
    def getCurrent(self,b):
        return self.currentData[b.id]
    def getCumulative(self,b):
        return self.cumulativeData[b.id]

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()],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,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='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ,PyRunner(iterPeriod=1,command="someCalculator.collect()",label='collect',dead=True)
 ,PyRunner(iterPeriod=1,command="someCalculator.pprint()",label='pprint',dead=True)
 ,PyRunner(iterPeriod=1,command="vtkExport()",label='vtkExp',dead=True)
 ]

someCalculator = SomeCalculator()

def triaxDone():
 triax.dead=True
 O.pause()

vtk = export.VTKExporter("heatData")
def vtkExport():
    vtk.exportSpheres(what=dict(current="someCalculator.getCurrent(b)",cumulative="someCalculator.getCumulative(b)"))

builtins.someCalculator = someCalculator

O.run(1000000,1)
O.step()

#### Applying normal stress
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  boundaryVel=copysign(min(0.01,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   print ('Normal stress =',currentSN,' | unbF=',unbF)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   print ('DONE! iter=',O.iter)
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  boundaryVel=copysign(min(0.1,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=10,label='servo')]+O.engines[4:]

O.run(1000000,1)

vtkExp.dead=False
collect.dead=False
pprint.dead=False
newton.damping=DAMPSHEAR
fixVel.dead = False

def fixVelocity(RATE_shear):
    O.bodies[0].state.vel[0] = RATE_shear
    topLayer_id = identify_topLayer()
    for i in topLayer_id:
        O.bodies[i].state.vel[0] = RATE_shear

O.run(100,1)
###

Thank you for taking your time!

Cheers,
Chien-Cheng

> I simplified my script
> O.run(1000000,1)

thanks, but is it possible to make it run shorter to demonstrate the error?

thanks
Jan

>is it possible to make it run shorter to demonstrate the error?

Yes, I also adjusted other parameters to make it run shorter to display the error.

I sent the modified script to you by email to avoid making the layout here too messy.
Hope you don't mind.

thanks!
Chien-Cheng

Thanks for more details

> self.cumulativeData = dict((a,0) for a in range(2,10))

here you tell that you only want to store bodies with ids in range(2,10)...

> vtk.exportSpheres(...)

... but exportSpheres exports all spheres by default, resulting in the error

You can (depending on your needs):

a) limit the export to the range:
vtk.exportSpheres(ids=list(range(2,10)),...)

b) return some artificial value for not-stored bodies:
class SomeCalculator:
   ...
   def getCurrent(self,b): # the same modification for getCumulative
         return self.currentData.get(b.id,0) # 0 is the "default" value for not-stored bodies, can be anything else

cheers
Jan

Dear Jan,

Sorry that I still have one more question...

If I would like to further use the cumulative data to calculate other parameters, such as temperature, how can I adjust the "singleton" class to print out the temperature?

thanks for your help again!
Chien-Cheng

to compute temperature, probably some loop inside collect method:
###
   def __init__(self):
      ...
      self.temperatureData = dict(self.cumulativeData) # one more copy
   ...
   def collect(self):
        self.currentData = ...
        self.computeTemperature() # here
   def computeTemperature(self): # compute and store actual temperature
      for k,cur in self.currentData.items():
         cum = self.cumulativeData[k]
         temerature = ... # some computation here
         self.temperatureData[k] = temperature
###

to print, just print it then :-)
###
    def pprint(self):
        ...
        print("temperatures")
        print("============")
        for key,valve in self.temperatureData.items():
            print(key,valve)
###

What I have created was just to show the idea of the approach, according to your specific needs, you can modify it to look completely different :-)

Also if the computation takes too long with Python, it is possible to create a C++ engine. But it is a nice feature that you can play easily with Python beforehand :-)

cheers
Jan

Hi Jan,

Thanks again for your reply and suggestion.

I am still learning programming with Python and am also a beginner in Yade. I am really thankful for these codes which save me plenty of time and allow me to learn faster. I will try to modify it with my own needs after getting more familiar with the details.

I tried your code but it turned out the cumulative and temperature cannot be printed out...

###
class SomeCalculator:
    def __init__(self):
        self.cumulativeHeatData = dict((a,0) for a in range(2,10))
        self.currentHeatData = dict(self.cumulativeHeatData)
        self.temperatureData = dict(self.cumulativeHeatData)
    def collect(self):
        for key,value in self.currentHeatData.items():
           self.cumulativeHeatData[key] += value
        self.currentHeatData = dict((a,self.getHeatFromOneParticle(a)) for a in range(2,10))
        self.computeTemperature()
    def getHeatFromOneParticle(self,a):
        b = O.bodies[a]
        totalHeat = 0
        for i in b.intrs():
            O.interactions.erase(i.id1,0)
            O.interactions.erase(i.id1,1)
            O.interactions.erase(0,i.id2)
            O.interactions.erase(1,i.id2)
            if not i.isReal: continue
            penetrationDepth = i.geom.penetrationDepth
            radius1 = i.geom.refR1
            radius2 = i.geom.refR2
            effectiveRadius = (radius1*radius2)/(radius1+radius2)
            contactRadius = math.sqrt(penetrationDepth*effectiveRadius)
            contactArea = np.pi*(contactRadius**2)
            relativeVelocity = i.geom.shearInc.norm()/O.dt
            generatedHeat = (i.phys.shearForce.norm()/contactArea)/(i.phys.normalForce.norm()/contactArea) * i.phys.normalForce.norm() * relativeVelocity
            totalHeat = totalHeat + generatedHeat
        return totalHeat
    def computeTemperature(self):
        heatCapacity = 900 # [J/kg*K]
        density = 2500 # [kg/m3]
        for key,cur in self.currentHeatData.items():
            cum = self.cumulativeHeatData[key]
            for a in range(2,10):
                b = O.bodies[a]
                temperature = 0.5 * cum * O.dt / heatCapacity / density * (4/3 * np.pi * b.shape.radius**3)
                self.temperatureData[key] = temperature
    def pprint(self):
        print()
        print("current heat")
        print("============")
        for key,value in self.currentHeatData.items():
            print(key,value)
        print("cumulative heat")
        print("==========")
        for key,value in self.cumulativeHeatData.items():
            print(key,value)
        print("temperatures")
        print("============")
        for key,valve in self.temperatureData.items():
            print(key,valve)
###

Is there something wrong with my code?

Cheers,
Chien-Cheng

Hello,

I have tried your code (put SomeClass in the script you have sent me) and the temperature was printed in terminal without problem..

> vtk.exportSpheres(what=dict(current="someCalculator.getCurrent(b)",...

it was not saved to vtk, because of the error
AttributeError: 'SomeCalculator' object has no attribute 'getCurrent'

If it behaves differently on your side, please provide more info

cheers
Jan

Hi Jan,

My mistake... the iterPeriod value I input for "someCalculator.collect()" and "someCalculator.pprint()" is different.
That explains the reason why I saw the zero value in cumulative data and temperature.

Cheers,
Chien-Cheng

Thanks Jan Stránský, that solved my question.