Outputs of the Triaxial code by Bruno Chareyre

Asked by ehsan benabbas on 2020-01-09

Hi everyone,

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

I use the Triaxial code by Bruno Chareyre [1] to understand how it works.

My goal is to code a Triaxial test with 20000 particles in different sizes (sand) under different strain and stress paths and get all inter-particle forces, branch vectors, contact normals, and fabric tensor. (with the properties defined in the code)

Based on the code and its output which can be seen in the following, I have asked my questions at the end of this message.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

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

print ('============ DEFINING VARIABLES ============')

num_spheres=20000,
compFricDegree = 30,
key='_triax_base_',
unknownOk=True
)

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

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

print ('============ DEFINING FUNCTIONS ============')
def history():
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,
n = triax.porosity)

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

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

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

print ('============ DEFINING PACKING ============')
frictionlesswalls=aabbWalls([mn,mx],thickness=0,material='frictionlesswalls')
wallIds=O.bodies.append(frictionlesswalls)

sp=pack.SpherePack()
clumps=False
if clumps:
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeClumpCloud(mn,mx,[c1],periodic=False)
sp.toSimulation(material='spheres')
O.bodies.updateClumpProperties()
else:
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)

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

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

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

print ('============ DEFINING ENGINES ============')
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,
PyRunner(iterPeriod=20,command='history()',label='recorder'),
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
NewtonIntegrator(damping=damp,label="newton")
]
Gl1_Sphere.stripes=True
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.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 ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
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))
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]

######################################################

triax.internalCompaction=False
triax.goal1 = sigmaIso
triax.goal2 = rate
triax.goal3 = sigmaIso
newton.damping = 0.1
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 ~~~~~~~~~~~~~')
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',None,'ev')}
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$'}
plot.plot()
plot.saveDataTxt('results'+key)
plot.saveGnuplot('plotScript'+key)
print ('************** END **************')
O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now output:

ehsan@ehsan:~/Desktop$/home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 Q.py Welcome to Yade 2019-08-08.git-775ae74 Using python version: 3.6.9 (default, Nov 7 2019, 10:44:02) [GCC 8.3.0] TCP python prompt on localhost:9000, auth cookie scdyua' XMLRPC info provider on http://localhost:21000 Running script Q.py ************** START ************** ============ DEFINING VARIABLES ============ ============ DEFINING FUNCTIONS ============ ============ DEFINING MATERIALS ============ ============ DEFINING PACKING ============ ============ DEFINING TRIAXIAL TEST ============ ============ DEFINING ENGINES ============ The constructor with a shareWidget is deprecated, use the regular contructor instead. Number of elements: 20006 Box Volume: 6.90392657805617e-310 Box Volume calculated: 8000.0 ============ APPLYING CONFINING PRESSURE ============ ~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~ unbalanced force: 0.007284146090760937 mean stress engine: -1.3899805976256234 mean stress (Calculated): -2.537923714531715 porosity= 0.8343868539188172 void ratio= 5.038168005756043 ################## Isotropic phase is finished and saved successfully ################## ============ APPLYING DEVIATORIC LOADING ============ step= 1000 unbalanced force: 0.007284146090760937 sigma2: -0.0 q= 50000.0 axial deformation (%) 0.0 ~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~ ################## Deviatoric phase is finished and saved successfully ################## ============ RECORD AND PLOT DATA ============ /home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/plot.py:444: MatplotlibDeprecationWarning: The 'verts' kwarg was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use 'marker' instead. scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,s=scatterSize,color=line.get_color(),**scatterMarkerKw) ************** END ************** Analysis has been taken for 270.053 seconds or 4.500883333333333 minutes [[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]] In [1]: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Questions: 1- Is there any problem with getting the following message in the output right after "DEFINING ENGINES"? >> The constructor with a shareWidget is deprecated, use the regular contructor instead. 2- Why the volume box calculated by triax is near to 0? >> Box Volume: 6.90392657805617e-310 3- Why in the "APPLYING CONFINING PRESSURE" section the porosity is about 0.8 while I defined it as 0.35? >> porosity= 0.8343868539188172 Thank you, Ehsan Question information Language: English Edit question Status: Expired For: Yade Edit question Assignee: No assignee Edit question Last query: 2020-01-17 Last reply: 2020-02-01 This question was reopened  Robert Caulk (rcaulk) said on 2020-01-12: #1 Please run the uncommented example script without modifying it and tell us if everything works as intended. >> 1- Is there any problem with getting the following message in the output right after "DEFINING ENGINES"? >> The constructor with a shareWidget is deprecated, use the regular contructor instead. How did you install yade? Seems like a libqglviewer problem. >> 3- Why in the "APPLYING CONFINING PRESSURE" section the porosity is about 0.8 while I defined it as 0.35? >> porosity= 0.8343868539188172 Because you arent using the part of the example script that reaches the target porosity [1].  ehsan benabbas (ehsanben) said on 2020-01-12: #2 Thank you for your response. This is the uncommented script. from yade import pack nRead=readParamsFromTable( num_spheres=1000,# number of spheres compFricDegree = 30, # contact friction during the confining phase key='_triax_base_', # put you simulation's name here unknownOk=True ) from yade.params import table num_spheres=table.num_spheres# number of spheres key=table.key targetPorosity = 0.43 #the porosity we want for the packing compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process) finalFricDegree = 30 # contact friction during the deviatoric loading rate=-0.02 # loading rate (strain rate) damp=0.2 # damping coefficient stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below) young=5e6 # contact stiffness mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing ## create materials for spheres and plates O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres')) O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls')) ## create walls around the packing walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) ## use a SpherePack object to generate a random loose particles packing sp=pack.SpherePack() clumps=False #turn this true for the same example with clumps if clumps: ## approximate mean rad of the futur dense packing for latter use volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) mean_rad = pow(0.09*volume/num_spheres,0.3333) ## define a unique clump type (we could have many, see clumpCloud documentation) c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)]) ## generate positions and input them in the simulation sp.makeClumpCloud(mn,mx,[c1],periodic=False) sp.toSimulation(material='spheres') O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia else: sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp]) triax=TriaxialStressController( ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions. ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth) finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth) thickness = 0, ## switch stress/strain control using a bitmask. What is a bitmask, huh?! ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0. ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e. ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc. stressMask = 7, internalCompaction=True, # If true the confining pressure is generated by growing particles ) 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()] ), ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf) GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), newton ] #Display spheres with 2 colors for seeing rotations better Gl1_Sphere.stripes=0 if nRead==0: yade.qt.Controller(), yade.qt.View() triax.goal1=triax.goal2=triax.goal3=-10000  ehsan benabbas (ehsanben) said on 2020-01-12: #3 I run this, there is no problem. and this its output: ehsan@ehsan:~/Desktop$ /home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 net.py
Using python version: 3.6.9 (default, Nov 7 2019, 10:44:02)
[GCC 8.3.0]
TCP python prompt on localhost:9000, auth cookie ausdcy'
XMLRPC info provider on http://localhost:21000
Running script net.py
The constructor with a shareWidget is deprecated, use the regular contructor instead.
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

In [1]:

 ehsan benabbas (ehsanben) said on 2020-01-12: #4

I have installed yade source code based on the instruction. Do you think I should install libqglviewer again?

I used reaching the target porosity but it didn't converged . that is why i tried to run the code without this part and then add this part to the code.

 ehsan benabbas (ehsanben) said on 2020-01-12: #5

Please check the followings. It seemo the package has not been installed. right?

ehsan@ehsan:~$QOpenGLWidget --version QOpenGLWidget: command not found ehsan@ehsan:~$ dpkg -s QOpenGLWidget | grep 'Version'
dpkg-query: package 'qopenglwidget' is not installed and no information is available
Use dpkg --info (= dpkg-deb --info) to examine archive files,
and dpkg --contents (= dpkg-deb --contents) to list their contents.
ehsan@ehsan:~$which QOpenGLWidget ehsan@ehsan:~$

 Robert Caulk (rcaulk) said on 2020-01-13: #6

>>This is the uncommented script.

That is not the example tutorial script. The script you link to, ([1]) is the example triax tutorial script.

 ehsan benabbas (ehsanben) said on 2020-01-13: #7

The code I shared as the uncommented is exactly this one. I only separated the uncommented lines. We are talking about the same code.

 Robert Caulk (rcaulk) said on 2020-01-14: #8

>>We are talking about the same code.

Ok, the code you post in #2 is missing most of the functionality of [1]. Including reaching a specified porosity.

 ehsan benabbas (ehsanben) said on 2020-01-14: #9

I am using the "reaching porosity" part but it seems like never would be converged.

my target porosity is 0.4 with frictional angle 30.

It still running but the output is like the following for 1000 iterations

Using python version: 3.6.9 (default, Nov 7 2019, 10:44:02)
[GCC 8.3.0]
TCP python prompt on localhost:9000, auth cookie `cusdye'
XMLRPC info provider on http://localhost:21000
Running script Q.py
************** START **************
============ DEFINING VARIABLES ============
============ DEFINING FUNCTIONS ============
============ DEFINING MATERIALS ============
============ DEFINING PACKING ============
============ DEFINING TRIAXIAL TEST ============
============ DEFINING ENGINES ============
The constructor with a shareWidget is deprecated, use the regular contructor instead.
Number of elements: 20006
Box Volume: 6.92593712388204e-310
Box Volume calculated: 8000.0
============ APPLYING CONFINING PRESSURE ============
~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~
unbalanced force: 0.0007999999999999998 mean stress engine: -0.0012625089975651174 mean stress (Calculated): -0.0
porosity= 0.8740795009613841
void ratio= 6.941518717245006
################## Isotropic phase is finished and saved successfully ##################
============ REACHING TARGET POROSITY ============
Friction: 28.5 porosity: 0.8740795009613841
Friction: 27.075 porosity: 0.8621664930170699
Friction: 25.721249999999998 porosity: 0.8476063035032143
.
.
.
Friction: 0.01957026971200895 porosity: 0.482076166440355
Friction: 0.0185917562264085 porosity: 0.4777297594206758

So the frictional angle is going to be 0. Actually I am still wondering how to define a packing with 30 degree for frictional angle, "20mm*20mm*20mm" dimensions, initial porosity of 0.4, 2e4 number of particles, minimum radius=0.1, and maximum radius=0.3 because this is the sample I want to study.

 Robert Caulk (rcaulk) said on 2020-01-15: #10

Hello,

When I run the example triaxial tutorial script [1], it reaches the target porosity within 1 minute.

As requested in comment #1, please run the example triaxial tutorial script [1] without removing lines of code and without modifying parameters. Does it work as designed?

Cheers,

Robert

 ehsan benabbas (ehsanben) said on 2020-01-15: #11

Hello Robert an thank you for your asnwers

Yes I did run the [1] without making any changes and it works perfectly. it convergence fast and show the deformation really smooth. This is what I wish to get for my model. So, I just made the following changes:

num_spheres=20000
targetPorosity = 0.4
mn,mx=Vector3(0,0,0),Vector3(20,20,20)
poisson=0.5
sp.makeCloud(mn,mx,0.2,0.5,num_spheres,False, 0.95,seed=1)
triax.goal1=triax.goal2=triax.goal3=-0.3
if unb<stabilityThreshold and abs(-0.3-triax.meanStress)/0.3<0.001:
triax.goal1=-0.3
triax.goal3=-0.3

I only made these changes, each on its specific line. The run took about 2 days and it couldn't be converged in Confining section and also couldn't be converged to reach the porosity. the became almost 0 and porosity was fluctuating still.

 Jan Stránský (honzik) said on 2020-01-15: #12

> young=5e6
> triax.goal1=triax.goal2=triax.goal3=-0.3
> if unb<stabilityThreshold and abs(-0.3-triax.meanStress)/0.3<0.001:
>
> it couldn't be converged in Confining section and also couldn't be converged to reach the porosity

what is a typical value of triax.meanStress? is there any chance the second condition can be non-randomly reach?

> the became almost 0

the what?

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-01-15: #13

Thanks Jan for clarification.

> what is a typical value of triax.meanStress? is there any chance the second condition can be non-randomly reach?
There is no typical value of mean stress. It's an average value of stresses in X,Y, and Z directions. Yes maybe, I should check, however I think it did not converge.

>the what?
the frictional angle became 0. sorry for not typing that. this specific one is solved as Robert answered me in another question. Still the problem exists.

Yes I did run the [1] without making any changes and it works perfectly. it convergence fast and show the deformation really smooth. This is what I wish to get for my model. So, I just made the following changes:

num_spheres=20000
targetPorosity = 0.4
mn,mx=Vector3(0,0,0),Vector3(20,20,20)
poisson=0.5
sp.makeCloud(mn,mx,0.2,0.5,num_spheres,False, 0.95,seed=1)
triax.goal1=triax.goal2=triax.goal3=-0.3
if unb<stabilityThreshold and abs(-0.3-triax.meanStress)/0.3<0.001:
triax.goal1=-0.3
triax.goal3=-0.3

I only made these changes, each on its specific line. The run took about 2 days and it couldn't be converged in Confining section and also couldn't be converged to reach the porosity. the became almost 0 and porosity was fluctuating still.

Bests,
Ehsan

 Jan Stránský (honzik) said on 2020-01-15: #14

>> what is a typical value of triax.meanStress? is there any chance the second condition can be non-randomly reach?
> There is no typical value of mean stress. It's an average value of stresses in X,Y, and Z directions. Yes maybe, I should check, however I think it did not converge.

Sorry for not being specific enough. My doubts are about the value triax.goal1=-0.3.
Why did you choose this value? What is wrong with original 10000?
The typical value of triax.meanStress I meant w.r.t. the defined goal -0.3.
With young=5e6, I would expect that even random stress fluctuations would be much higher than 0.3. That's why I asked if the condition (-0.3-triax.meanStress)/0.3<0.001 could be fulfilled .

Jan

 ehsan benabbas (ehsanben) said on 2020-01-16: #15

Oh you're right. I didn't pay attention to "young". as my specimen dimensions are 20 mm. the confining stresses should be "300Kpa" which is equal to "0.3 Mpa= 0.3 N/mm2" . I chose Mpa to make "mm" consistent. So, in this case I think I should change the value of "young". I'm not sure what's the unit of "5e6 young" in the code. Could you please tell me its unit? I should make it consistent with Newton and mm. Even I should make the unit of density consistent. So in this case, since the density is equal to 2e3 kg/m3 , it should be change to 2e13 N/mm3 as well. right?

 ehsan benabbas (ehsanben) said on 2020-01-16: #16

For now, I change it to 300 Mpa.

 Jan Stránský (honzik) said on 2020-01-16: #17

Why not simply use basic units? then you are sure you are consistent and prevent confusion.

Yade has no units, it just computes numbers. Stress might be Pa, N/m2, MPa, N/mm2 as well as something like pound/inch2.
It was discussed several times before, maybe we can create a FAQ on this topic..

cheers
Jan

 Robert Caulk (rcaulk) said on 2020-01-16: #18

>> So, I just made the following changes

Thank you! This is the sort of concise information that allows us to see exactly what is "wrong" and how to fix it (thanks Jan ;) ).

>>maybe we can create a FAQ on this topic..

I am guessing it would be as futile as the how to ask ;).

 ehsan benabbas (ehsanben) said on 2020-01-17: #19

Jan and Robert, thanks for your feedback.

When I want to use basic units, it implicitly means that I should make them consistent in my mind to a specific unit. For example, make them consistent based on "mm". Exactly like what I did. Accordingly, I made changes for "m" and "N" as follows:

Sima iso = 300 Kpa = 3e5 N/m2
Porosity = 0.4
Nu = 0.5
Phi = 30
Density = 2000 Kg/m3
E = 5e7 pa

Which needs to make these changes in the code:

num_spheres=20000
targetPorosity = 0.4
young=5e7
mn,mx=Vector3(0,0,0),Vector3(0.02,0.02,0.02)
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
sp.makeCloud(mn,mx,0.0002,0.5,num_spheres,False, 0.95,seed=1)
triax.goal1=triax.goal2=triax.goal3=-300000
triax.goal1=-300000
triax.goal3=-300000

In this case, it takes about 2 mins to get the results. But plots are not conventional. Would you please make these changes and run the modified [1] and let me know your feedback?

 ehsan benabbas (ehsanben) said on 2020-01-17: #20

Sorry I still need an answer but seems wrongly chose the problem is solved

 Launchpad Janitor (janitor) said on 2020-02-01: #21

This question was expired because it remained in the 'Open' state without activity for the last 15 days.