The problem of plot

Asked by Zhicheng Gao

The first problem is that when I run my code, the following exceptions and errors occur:
Unhandled exception in thread started by <function liveUpdate at 0x7f9af8eb7350>
Traceback (most recent call last):
  File "/usr/lib/x86_64-linux-gnu/yade/py/yade/plot.py", line 508, in liveUpdate
    ax.relim() # recompute axes limits
  File "/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py", line 1938, in relim
    self._update_line_limits(line)
  File "/usr/lib/python2.7/dist-packages/matplotlib/axes/_base.py", line 1801, in _update_line_limits
    path = line.get_path()
  File "/usr/lib/python2.7/dist-packages/matplotlib/lines.py", line 957, in get_path
    self.recache()
  File "/usr/lib/python2.7/dist-packages/matplotlib/lines.py", line 667, in recache
    self._xy = np.column_stack(np.broadcast_arrays(x, y)).astype(float)
  File "/usr/lib/python2.7/dist-packages/numpy/lib/stride_tricks.py", line 249, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/usr/lib/python2.7/dist-packages/numpy/lib/stride_tricks.py", line 184, in _broadcast_shape
    b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape

The second problem: I use the plot.resetData() and plot.live= True to clear the previous data and update the plot, but it doesn't work, the plot isn't updated.

The last problem is that I use the O.cell.trsf=Matrix3.Identity to set the current cell configuration to be the reference one, but when I use print(O.cell.trsf, triax.strain) to print the current strain, the strain matrix is not the identity matrix.

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Janek Kozicki
Solved:
Last query:
Last reply:
Revision history for this message
Zhicheng Gao (zhichenggao) said :
#1

this is my code:
##______________ First section, generate sample_________

from __future__ import print_function
from yade import pack, qt, plot
from math import *

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=100,
        targetPorosity= .387,
        confiningPressure=-100000,
        ## material parameters
        compFricDegree=15,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2600,
        alphaKr=7.5,
        alphaKtw=0,
 competaRoll=.22,
        finaletaRoll=.22,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

O.periodic=True
O.cell.hSize=Matrix3(.001,0,0, 0,.001,0, 0,0,.001)
# create materials for spheres
#shear strength is the sum of friction and adhesion, so the momentRotationLaw=True
O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres'))

# generate particles packing
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.001,.001,.001),psdSizes=[0.00008,0.000125,0.0001592,0.0002003,0.0003153,0.000399,0.000502,0.0005743],psdCumm=[0.0,0.00628,0.0565,0.198,0.721,0.915,0.991,1.0],num=num_spheres,seed=1)
sp.toSimulation(material='spheres')

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D()],
            [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)],
            [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)],
 ),
        PeriodicFlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        PeriTriaxController(label='triax',
            # specify target values and whether they are strains or stresses
            goal=(confiningPressure,confiningPressure,confiningPressure), stressMask=7,
            # type of servo-control, the strain rate isn't determined, it shloud check the unbalanced force
            dynCell=True,maxStrainRate=(10,10,10),
            # wait until the unbalanced force goes below this value
            maxUnbalanced=stabilityThreshold,relStressTol=1e-3,
            doneHook='compactionFinished()'
            ),
        NewtonIntegrator(damping=0)
]
qt.View()

import sys
def compactionFinished():
        #check the current porosity
        # if the current porosity is lager than target Porosity and comFricDegree is lager than 10,
        # then we decrease friction value and apply it to all the bodies and contacts,
        # else we decrease rolling friction value.
 global compFricDegree, competaRoll
        if porosity()>targetPorosity and compFricDegree>5:
            # we decrease friction value and apply it to all the bodies and contacts
            compFricDegree = 0.95*compFricDegree
            setContactFriction(radians(compFricDegree))
            print('Friction:', compFricDegree,'porosity:', porosity())
            # python syntax, make each step printout
            sys.stdout.flush()
        elif porosity()>targetPorosity:
            # we decrease rolling fiction value and apply it to all the bodies and contacts
     competaRoll=0.95*competaRoll
            for b in O.bodies:
                b.mat.etaRoll=competaRoll
            for i in O.interactions:
                i.phys.etaRoll=competaRoll
            print('Rollingfriction:', b.mat.etaRoll, 'porosity:', porosity())
            sys.stdout.flush()
        else:
     # after sample preparation, save the state
            O.save('compactedState'+filename+'.yade.gz')
            print('Compacted state saved', 'porosity', porosity())
            # next time, called python command
            triax.doneHook=''
            O.pause()

# enable energy tracking in the code
O.trackEnergy=True

# define function to record history
def history():
    plot.addData(unbalanced=unbalancedForce(),i=O.iter,exx=-triax.strain[0],
            eyy=-triax.strain[1], ezz=-triax.strain[2],
            sxx=-triax.stress[0],syy=-triax.stress[1],szz=-triax.stress[2],
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            porosity=porosity(),Etot=O.energy.total(),**O.energy# save all available energy data
            )
O.engines=O.engines+[PyRunner(command='history()', iterPeriod=20)]
# define what to plot
plot.plots={'i':('unbalanced','porosity'),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),' i ':('Etot')}
# show the plot
plot.plot()
plot.live=True

O.run()
O.wait()

plot.saveDataTxt('compactedState'+filename)
##__________________________________second section, deviatoric loading__________________________
# change the contact parameters to the final calibration value
setContactFriction(radians(finalFricDegree))
for b in O.bodies:
    b.mat.etaRoll=finaletaRoll
for i in O.interactions:
    i.phys.etaRoll=finaletaRoll
print(O.cell.hSize,O.cell.trsf)
# set the current cell configuration to be the reference one
O.cell.trsf=Matrix3.Identity
print(O.cell.trsf, triax.strain)
# change control type: keep constant confinement in x,y, 30% compression in z
triax.goal=(confiningPressure,confiningPressure,-.3)
triax.stressMask=3
# allow faster deformation along x,y to better maintain stresses
triax.maxStrainRate=(1.,1.,.01)
# call triaxFinished instead of compactionFinished
triax.doneHook='triaxFinished()'

def triaxFinished():
        O.save('loadedState'+filename+'.yade.gz')
 print('Finished')
 O.pause()

# Reset all plot data; keep plots and labels intact
plot.resetData()
O.run(1000,True)
plot.saveDataTxt('loadedState'+filename)

Revision history for this message
Janek Kozicki (cosurgi) said :
#2
Revision history for this message
Janek Kozicki (cosurgi) said :
#3

"ValueError: shape mismatch" this happens because a plot uses two lists of points: arguments (say "x") and values (say "y"). Plotting happens at the same time when calculations are pushing more data points to the plot. Plotting fails when the amount of "x" differs from the amount of "y" which happens quite often if calcuations are adding "x" then "y" in a separate thread. And when their sizes differ then it can't be plotted. This function setLiveForceAlwaysUpdate makes sure that calculations are paused during plot refresh, so the problem won't happen.

This bug was reported in https://gitlab.com/yade-dev/trunk/-/issues/110 and fixed in https://gitlab.com/yade-dev/trunk/-/merge_requests/570/diffs

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4

Thanks, Janek. I still have a problem, I use the O.cell.trsf=Matrix3.Identity to set the current cell configuration to be the reference one, but when I use print(O.cell.trsf, triax.strain) to print the current strain, the strain is not equal to zero.

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#5

Dear Janek,
  According to your suggestion, I add plot.setLiveForceAlwaysUpdate(True) to my code, but the following error occurred:
Traceback (most recent call last):
  File "/usr/bin/yade", line 182, in runScript
    execfile(script,globals())
  File "plottest.py", line 136, in <module>
    plot.setLiveForceAlwaysUpdate(True)
AttributeError: 'module' object has no attribute 'setLiveForceAlwaysUpdate'

Revision history for this message
Best Janek Kozicki (cosurgi) said :
#6

AttributeError: 'module' object has no attribute 'setLiveForceAlwaysUpdate'

You need to upgrade your version to 2021.01a, this was fixed only 5 months ago.

The other one is a different question, you need to open a new question for that.

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#7

Thanks Janek Kozicki, that solved my question.

Revision history for this message
doclarens (doclarens) said :
#8

hanks, Janek. I still have a problem, I use the AttributeError: 'module' object has no attribute 'setLiveForceAlwaysUpdate' to set the current cell configuration to be the reference one, but when I use print File "plottest.py", line 136, in <module> to print the current strain, the strain is not equal to zero. https://dentnis.com