PeriTriaxController Load

Asked by Johannes Welsch on 2021-03-02

Hi I'm having a little problem with the PeriTriaxController and saving and loading of simulations. Maybe one of you can help me.
I provided a MWE below.

My Problems are:
 - When I first run the simulation everything works fine, but simulation doesn't stop when the strain -0.4 is reached, so I think the doneHook isn't triggered, but I don't know why
 - When I rund the simulation for the second time, there is already the file 'three.yade.gz' containing the information of the first run. The problem here is that the shearing of the sample doesn't start. When looking at the output, it goes to phase 2, but there nothing happens.

Thank you in advance for your help.
Johannes

# encoding: utf-8
from yade import pack, qt, plot

sigmaIso=-25000

Filename = 'three.yade.gz'
savedState = os.path.exists(Filename)

if not savedState:
    O.periodic=True
    O.cell.setBox(5,5,10)

    compFricDegree = 35
sphere=O.materials.append(FrictMat(young=1e7,poisson=0.5,density=2650,frictionAngle=radians(compFricDegree),label='sphere'))
    psdSizes,psdCumm = [0.99,1.01],[0,1.0]

    sp=pack.SpherePack()
    sp.makeCloud(Vector3(0,0,0),O.cell.refSize,psdSizes=psdSizes,psdCumm=psdCumm,num=1000,seed=1,distributeMass=True,periodic=True)
    sp.toSimulation()
else:
    O.load(Filename)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],
                    [Ip2_FrictMat_FrictMat_FrictPhys()],
                    [Law2_ScGeom_FrictPhys_CundallStrack()]
                    ),
    PeriTriaxController( label='triax',
                            goal=(sigmaIso,sigmaIso,sigmaIso),
                            stressMask=7,
                            dynCell=True,
                            maxStrainRate=(1000,1000,1000),
                            maxUnbalanced=1e-3,
                            relStressTol=1e-3,
                            doneHook='triaxFinished()'
                        ),
    NewtonIntegrator(damping=.1),
]

O.dt=.5*PWaveTimeStep()

def triaxFinished():

    print 'Making decision'

    if not savedState:
        print 'Going to Phase 1'
        triax.doneHook = 'Phase_1()'
    else:
        print 'Going to Phase 2'
        triax.doneHook = 'Phase_2()'

def Phase_2():
    print 'Phase 2'
    O.cell.setBox(O.cell.size[0],O.cell.size[1],O.cell.size[2])
    triax.goal = (sigmaIso,sigmaIso,-0.4)
    triax.stressMask = 3
    triax.maxStrainRate = (0,0,0.05)
    triax.maxUnbalanced = 0.1
    triax.doneHook='Stop()'

def Phase_1():
    print 'Phase 1'
    O.save(Filename)
    print 'Saved'
    triax.doneHook = 'Phase_2()'

def Stop():
    O.pause()

def UnbalancedForce():
    print unbalancedForce(), triax.strain[2]

O.engines += [PyRunner(command='UnbalancedForce()',iterPeriod=1000)]

O.run()
O.SaveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Johannes Welsch
Solved:
2021-03-02
Last query:
2021-03-02
Last reply:
2021-03-02
Jan Stránský (honzik) said : #1

Hello,

> When I first run the simulation everything works fine, but simulation doesn't stop when the strain -0.4 is reached, so I think the doneHook isn't triggered, but I don't know why
> print 'Making decision'

is this printed?
If no, the goal is defined by stress, so strain -0.4 is not relevant (yet)

> triax.maxStrainRate = (0,0,0.05)

why zeros for x and y?

> O.engines += [PyRunner(command='UnbalancedForce()',iterPeriod=1000)]

Since there is no running between O.engines definition and this line, put the PyRunner directly to O.engines definition

cheers
Jan

Hello Jan,
thanks for your quick response

> O.engines += [PyRunner(command='UnbalancedForce()',iterPeriod=1000)]
I put it to PyRunner directly

> triax.maxStrainRate = (0,0,0.05)
it was zero for no reason, I changed it to 0.05 as in z-direction

> print 'Making decision'
gives the following output for first simulation (no savedState)

>Making decision
>Going to Phase 1
>Phase 1
>Saved
>Phase 2

and gives the output for second simulation (savedState is True)
>Making decision
>Going to Phase 2
>Making decision
>Going to Phase 2
... and so on for every iteration.
I really don't know why there seems to be nothin happening in Phase 2.

Hello,

I think I found the solution.
I added the definition of O.engines to the if clause of saved State as it is also saved in the yade.gz file. Now everything works properly.

if not savedState:
    O.engines=[
        ...
    ]

I also tried to output len(O.engines) in UnbalancedForce() at every 1000 iterations, but i wasn't showing me that I loaded the engines twice. It was still showing 6 engines. When I added O.engines = [..] in the if clause everything works.

Thank you very much for your help Jan Stransky

Here the working MWE:

# encoding: utf-8
from yade import pack, qt, plot

sigmaIso=-25000

Filename = 'three.yade.gz'
savedState = os.path.exists(Filename)

if not savedState:
    O.periodic=True
    O.cell.setBox(5,5,10)

    compFricDegree = 35

    sphere=O.materials.append(FrictMat(young=1e7,poisson=0.5,density=2650,frictionAngle=radians(compFricDegree),label='sphere'))
    psdSizes,psdCumm = [0.99,1.01],[0,1.0]

    sp=pack.SpherePack()
    sp.makeCloud(Vector3(0,0,0),O.cell.refSize,psdSizes=psdSizes,psdCumm=psdCumm,num=500,seed=1,distributeMass=True,periodic=True)
    sp.toSimulation()
else:
    O.load(Filename)

if not savedState:
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],
                        [Ip2_FrictMat_FrictMat_FrictPhys()],
                        [Law2_ScGeom_FrictPhys_CundallStrack()]
                        ),
        PeriTriaxController( label='triax',
                                goal=(sigmaIso,sigmaIso,sigmaIso),
                                stressMask=7,
                                dynCell=True,
                                maxStrainRate=(1000,1000,1000),
                                maxUnbalanced=1e-3,
                                relStressTol=1e-3,
                                doneHook='triaxFinished()'
                            ),
        PyRunner(command='UnbalancedForce()',iterPeriod=100),
        NewtonIntegrator(damping=.1),
    ]

O.dt=.5*PWaveTimeStep()

def triaxFinished():

    print 'Making decision'

    if not savedState:
        print 'Going to Phase 1'
        triax.doneHook = 'Phase_1()'
    else:
        print 'Going to Phase 2'
        triax.doneHook = 'Phase_2()'

def Phase_2():
    print 'Phase 2'
    O.cell.setBox(O.cell.size[0],O.cell.size[1],O.cell.size[2])
    triax.goal = (sigmaIso,sigmaIso,-0.2)
    triax.stressMask = 3
    triax.StrainRate = (0,0,-0.01)
    triax.maxStrainRate = (0.05,0.05,0.05)
    triax.maxUnbalanced = 0.01

    print unbalancedForce(), triax.strain[2]
    print triax.stress[2]
    print triax.stressMask
    print triax.strainRate
    #print triax.doneHook

    triax.doneHook='Stop()'

def Phase_1():
    print 'Phase 1'
    O.save(Filename)
    print 'Saved'
    triax.doneHook = 'Phase_2()'

def Stop():
    print 'Stop'
    O.pause()

def UnbalancedForce():
    print unbalancedForce(), triax.strain[2], len(O.engines)

O.run()

Jan Stránský (honzik) said : #5

Nice catch :-) one has to be careful when loading saved state to the same script it was saved from..
Jan

Note: I found that:

1/ the example for save/reload on the wiki page [1] was obsolete (for some reason the github.com link stopped working), it is fixed now [2].

2/ in the example the engines are re-defined. From this discussion I suspect it gives room to more issues than just keeping engines definition udner 'savedState' guard. Agreed?

B

[1] https://www.yade-dem.org/wiki/Howtoask (5th bullet point of 2nd section)
[2] https://gitlab.com/yade-dev/trunk/-/blob/3e9a209234b7f23241d5f4bdef1b586056e97582/examples/simple-scene/save-then-reload.py