yade-users team mailing list archive

[ehsan benabbas <question687726@xxxxxxxxxxxxxxxxxxxxx>]

New question #687726 on Yade:
https://answers.launchpad.net/yade/+question/687726

Hello all and happy new year!

I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74

I want to code a 3D-triaxial test problem which applies a rate of strain for the deviatoric part. I run my code which takes a lot of time, I left it to be run for 3 days, and get this error:

Obviously it cannot be converged in the deviatoric section!

********** ERROR ************

axial deformation (%) -2575.3325503377832
~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~
step= 882430 unbalanced force: 2.6340108400409105e-06 sigma2: -3927217.501133026 q= -3877217.501133026
axial deformation (%) -2575.3325503377832
~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~
step= 882500 unbalanced force: 2.4289637626032886e-06 sigma2: -3927217.5011330247 q= -3877217.5011330247
axial deformation (%) -2575.3325503377832
~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~
Traceback (most recent call last):
  File "/home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74", line 336, in runScript
    execfile(script,globals())
  File "/usr/lib/python3/dist-packages/past/builtins/misc.py", line 82, in execfile
    exec_(code, myglobals, mylocals)
  File "Q.py", line 154, in <module>
    O.run(1000, True)
RuntimeError: boost::thread_resource_error: Resource temporarily unavailable
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

^[[55;1R---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
~/yade/install/bin/yade-2019-08-08.git-775ae74 in <module>
    402 import yade.qt
    403 qapp=QApplication(sys.argv)
--> 404 userSession(gui=gui,qapp=qapp)
    405
    406 if __name__ == "__main__":

~/yade/install/bin/yade-2019-08-08.git-775ae74 in userSession(gui, qapp)
    351 if yade.runtime.opts.exitAfter: sys.exit(0)
    352 ## start interpreter
--> 353 if not opts.mpi_mode: ipshell()
    354
    355

~/.local/lib/python3.6/site-packages/IPython/terminal/embed.py in __call__(self, header, local_ns, module, dummy, stack_depth, global_ns, compile_flags, **kw)
    227 # our call and get the original caller's namespaces.
    228 self.mainloop(local_ns, module, stack_depth=stack_depth,
--> 229 global_ns=global_ns, compile_flags=compile_flags)
    230
    231 self.banner2 = self.old_banner2

~/.local/lib/python3.6/site-packages/IPython/terminal/embed.py in mainloop(self, local_ns, module, stack_depth, display_banner, global_ns, compile_flags)
    322
    323 with self.builtin_trap, self.display_trap:
--> 324 self.interact()
    325
    326 # now, purge out the local namespace of IPython's hidden variables.

~/.local/lib/python3.6/site-packages/IPython/terminal/interactiveshell.py in interact(self, display_banner)
    479
    480 try:
--> 481 code = self.prompt_for_code()
    482 except EOFError:
    483 if (not self.confirm_exit) \

~/.local/lib/python3.6/site-packages/IPython/terminal/interactiveshell.py in prompt_for_code(self)
    408 default=default,
    409 # pre_run=self.pre_prompt,# reset_current_buffer=True,
--> 410 **self._extra_prompt_options())
    411 return text
    412

~/.local/lib/python3.6/site-packages/prompt_toolkit/shortcuts/prompt.py in prompt(self, message, editing_mode, refresh_interval, vi_mode, lexer, completer, complete_in_thread, is_password, key_bindings, bottom_toolbar, style, color_depth, include_default_pygments_style, style_transformation, swap_light_and_dark_colors, rprompt, multiline, prompt_continuation, wrap_lines, enable_history_search, search_ignore_case, complete_while_typing, validate_while_typing, complete_style, auto_suggest, validator, clipboard, mouse_support, input_processors, reserve_space_for_menu, enable_system_prompt, enable_suspend, enable_open_in_editor, tempfile_suffix, inputhook, async_, default, accept_default, pre_run)
    736 return ensure_future(run_async())
    737 else:
--> 738 return run_sync()
    739
    740 @property

~/.local/lib/python3.6/site-packages/prompt_toolkit/shortcuts/prompt.py in run_sync()
    725 with self._auto_refresh_context():
    726 self.default_buffer.reset(Document(default))
--> 727 return self.app.run(inputhook=self.inputhook, pre_run=pre_run2)
    728
    729 def run_async():

~/.local/lib/python3.6/site-packages/prompt_toolkit/application/application.py in run(self, pre_run, set_exception_handler, inputhook)
    707 loop.set_exception_handler(handle_exception)
    708 try:
--> 709 return run()
    710 finally:
    711 loop.set_exception_handler(previous_exc_handler)

~/.local/lib/python3.6/site-packages/prompt_toolkit/application/application.py in run()
    680 def run():
    681 f = self.run_async(pre_run=pre_run)
--> 682 run_until_complete(f, inputhook=inputhook)
    683 return f.result()
    684

~/.local/lib/python3.6/site-packages/prompt_toolkit/eventloop/defaults.py in run_until_complete(future, inputhook)
    121 Return the Future's result, or raise its exception.
    122 """
--> 123 return get_event_loop().run_until_complete(future, inputhook=inputhook)

~/.local/lib/python3.6/site-packages/prompt_toolkit/eventloop/posix.py in run_until_complete(self, future, inputhook)
     64
     65 while not future.done():
---> 66 self._run_once(inputhook)
     67
     68 # Run one last time, to flush the pending `_calls_from_executor`s.

~/.local/lib/python3.6/site-packages/prompt_toolkit/eventloop/posix.py in _run_once(self, inputhook)
     83 " True when there is input ready. The inputhook should return control. "
     84 return self._ready_for_reading(None if wait else 0) != []
---> 85 self._inputhook_context.call_inputhook(ready, inputhook)
     86
     87 # Wait until input is ready.

~/.local/lib/python3.6/site-packages/prompt_toolkit/eventloop/inputhook.py in call_inputhook(self, input_is_ready_func, inputhook)
     76 os.write(self._w, b'x')
     77
---> 78 threading.Thread(target=thread).start()
     79
     80 # Call inputhook.

/usr/lib/python3.6/threading.py in start(self)
    844 _limbo[self] = self
    845 try:
--> 846 _start_new_thread(self._bootstrap, ())
    847 except Exception:
    848 with _active_limbo_lock:

RuntimeError: can't start new thread
Fatal Python error: This thread state must be current when releasing

Thread 0x00007f6c09bd3740 (most recent call first):

Thread 0x00007f6c09bd3740 (most recent call first):

Current thread 0x00007f6ba0d31700 (most recent call first):
  File "/usr/lib/python3.6/bz2.py", line 245 in write
  File "/home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/plot.py", line 663 in saveDataTxt
  File "Q.py", line 84 in checkUnbalanced
  File "<string>", line 1 in <module>

Thread 0x00007f6bdffff700 (most recent call first):
  File "/usr/lib/python3.6/selectors.py", line 376 in select
  File "/usr/lib/python3.6/socketserver.py", line 236 in serve_forever
  File "/home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/remote.py", line 142 in run

Thread 0x00007f6be4f9e700 (most recent call first):
  File "/usr/lib/python3.6/selectors.py", line 376 in select
  File "/usr/lib/python3.6/socketserver.py", line 236 in serve_forever
  File "/home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/remote.py", line 142 in run

Thread 0x00007f6c09bd3740 (most recent call first):
Aborted (core dumped)
ehsan@ehsan:~/Desktop$

******************************************************************************
******************************************************************************

and my code is as follows:

******************************* MY CODE *************************************

############################################################################################################################
######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS #########
############################################################################################################################

print ('************** START **************')
import numpy as np
import math
from yade import pack, plot, qt, export, utils
from datetime import datetime

######################################################
######### DEFINING VARIABLES #########

print ('============ DEFINING VARIABLES ============')
nRead=readParamsFromTable(
 num_spheres=20000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)

from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.43
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=15e6
poisson=0.4
mn,mx=Vector3(0,0,0),Vector3(20,20,20)
sigmaIso=-50e3

######################################################
######### DEFINING MATERIALS #########

print ('============ DEFINING MATERIALS ============')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionlesswalls'))

####################################################
######### DEFINING PACKING #########

print ('============ DEFINING PACKING ============')
walls=aabbWalls([mn,mx],thickness=0,material='frictionlesswalls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)
clumps=False
if clumps:
    volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
    mean_rad = pow(0.09*volume/num_spheres,0.3333)
    c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
    sp.makeClumpCloud(mn,mx,[c1],periodic=False)
    sp.toSimulation(material='spheres')
    O.bodies.updateClumpProperties()
else:
    sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)
    O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

##########################################################
######### DEFINING TRIAXIAL TEST #########

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False
)

######################################################
######### DEFINING FUNCTIONS #########

print ('============ DEFINING FUNCTIONS ============')
def checkUnbalanced():
 if unbalancedForce()<.05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
def history():
 plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
   ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
   i=O.iter,unbalanced=unbalancedForce()
)

####################################################
######### DEFINING ENGINES #########

print ('============ DEFINING ENGINES ============')
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton,
 PyRunner(iterPeriod=20,command='history()',label='recorder'),
 PyRunner(command='checkUnbalanced()',realPeriod=2)
]
Gl1_Sphere.stripes=True
if nRead==0: yade.qt.Controller(), yade.qt.View()
print ('Number of elements: ', len(O.bodies))
print ('Box Volume: ', triax.boxVolume)
print ('Box Volume calculated: ', volume)

###############################################################
######### APPLYING CONFINING PRESSURE #########

print ('============ APPLYING CONFINING PRESSURE ============')
triax.internalCompaction=False
setContactFriction=radians(compFricDegree)
triax.stressmask=7
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
    print ('unbalanced force:',unb,' mean stress engine: ',triax.meanStress,' mean stress (Calculated): ',meanS)
    print ('porosity=',triax.porosity)
    print ('void ratio=',triax.porosity/(1-triax.porosity))
    print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
    if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
        break
O.save('confinedPhase'+key+'.xml')
print ('################## Isotropic phase is finished and saved successfully ##################')
e22Check=triax.strain[1]
######################################################
######### DEVIATORIC LOADING #########

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction=radians(compFricDegree)
triax.stressMask = 5
triax.goal2=rate
triax.goal1=sigmaIso
triax.goal3=sigmaIso
while 1:
    O.run(1000, True)
    unb=unbalancedForce()
    axialS=triax.stress(triax.wall_top_id)[1]
    print ('step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-sigmaIso)
    print ('axial deformation (%)', (triax.strain[1]-e22Check)*100)
    print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
    if abs((triax.strain[1]-e22Check)-0.01)<=0.001:
        O.save('firstpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.02)<=0.001:
        O.save('secondpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.03)<=0.001:
        O.save('thirdpoint.xml')
    if abs((triax.strain[1]-e22Check)-0.05)<=0.001:
        O.save('fourthpoint.xml')
    if triax.strain[1]-e22Check>=0.3 :
        break
O.save('final.xml')
print ('################## Deviatoric phase is finished and saved successfully ##################')

########################################################
######### RECORD AND PLOT DATA #########

print ('============ RECORD AND PLOT DATA ============')
O.run(5000,True)
plot.plots={'e22':('s11','s22','s33'),'e22':('ev')}
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e11':'$\epsilon_{11}$' , 'e22':'$\epsilon_{22}$' , 'e33':'$\epsilon_{33}$' , 'ev':'$\epsilon_{V}$'}
plot.plot()
plot.saveDataTxt('results'+key)
plot.saveGnuplot('plotScript'+key)
print ('************** END **************')

******************************************************************************

So I have following questions:
1 - Why it takes so long to be run?
2- Is there any way to run this code in a shorter time?
3- Why it cannot be converged?
4- What is the meaning of this error and what do I do to debug it?

Thank you very much for you kind helps.
Bests,

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.