Optimise material parameter

Asked by Mithushan Soundaranathan on 2021-06-03

Hi,

I am optimising k1 for the luding material, by using o.reset in a for loop. The code is not going through the loop, the code just goes to the last iteration. For example, I am going iterate through 100 iteration and I set in a loop, "for i in range (0,n_iteration,1):". The for loop runs by it self initially and even before simulation starts, and i set to 99 when the simulation starts.

How can I run my simulation in the loop, for n number of iteration. And changing the parameter for each iteration.

Best,
Mithu

Here is my code

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d

o = Omega()

# Physical parameters
fr = 0.54
rho = 1550
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 50000
kp = 5.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.005
Cyl_height=0.045
cross_area=math.pi*(Tab_rad**2)

loading_exp = pd.read_csv (r'Gao_data_loading.csv')
loading_exp = pd.DataFrame(loading_exp, columns= ['porosity','compression_pressure'])
f1_loading = interp1d(loading_exp.porosity, loading_exp.compression_pressure, kind= 'cubic')
porosity_loading=np.linspace(0.65, 0.23, num=60, endpoint=True)
interpol_loading_exp=f1_loading(porosity_loading)

unloading_exp = pd.read_csv (r'Gao_data_loading.csv')
unloading_exp = pd.DataFrame(unloading_exp, columns= ['porosity','compression_pressure'])
f1_unloading = interp1d(unloading_exp.porosity, unloading_exp.compression_pressure, kind= 'cubic')
porosity_unloading=np.linspace(0.23, 0.32, num=60, endpoint=True)
interpol_unloading_exp=f1_unloading(porosity_unloading)

Comp_press= max(loading_exp.compression_pressure)
Comp_force=Comp_press*cross_area

n_iter=80
RMSE_save=[]
#*************************************
for i in range (0,n_iter,1):
    o.reset()
    compression_profile_save=[]
    # Add material
    mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

    # Spheres for compression

    sp=pack.SpherePack()
    sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-18*Diameter),(3.0*Diameter,3.0*Diameter,18.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=527)
    sp.toSimulation(material=mat1)
    walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

    # Add engines
    o.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                            Bo1_Wall_Aabb(),
                            Bo1_Facet_Aabb()
                            ]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
        Ig2_Facet_Sphere_ScGeom(),
        Ig2_Wall_Sphere_ScGeom()],
        [Ip2_LudingMat_LudingMat_LudingPhys()],
        [Law2_ScGeom_LudingPhys_Basic()]
    ),
    NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
    PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
    #DeformControl(label="DefControl")
    ]

    def checkForce():
        if O.iter < 10000:
            return
        highSphere = 0.0
        for b in O.bodies:
            if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
                highSphere = b.state.pos[2]
            else:
                pass

        O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
        global plate
        plate = O.bodies[-1]
        plate.state.vel = (0, 0, -10)
        O.engines = O.engines + [PyRunner(command='checkdata()', iterPeriod=1500)]
        fCheck.command = 'unloadPlate()'

    def unloadPlate():
        if abs(O.forces.f(plate.id)[2]) > Comp_force:
            plate.state.vel *= -1
            fCheck.command = 'stopUnloading()'

    def stopUnloading():
        if abs(O.forces.f(plate.id)[2]) == 0:
            o.pause()
            compression_profile=pd.DataFrame(compression_profile_save,columns=['porosity','Comp_pressure'])
            idx_max=np.argmax(compression_profile.Comp_pressure)
            loading_sim=compression_profile[:(idx_max+1)]
            f1_loading_sim = interp1d(loading_sim.porosity, loading_sim.Comp_pressure, kind= 'cubic')
            interpol_loading_sim=f1_loading_sim(porosity_loading)
            RMSE_loading=math.sqrt((sum((interpol_loading_exp-interpol_loading_sim)**2))/len(porosity_loading))
            unloading_sim=compression_profile[(idx_max+2):-1]
            f1_unloading_sim = interp1d(unloading_sim.porosity, unloading_sim.Comp_pressure, kind= 'cubic')
            interpol_unloading_sim=f1_unloading_sim(porosity_unloading)
            RMSE_unloading=math.sqrt((sum((interpol_unloading_exp-interpol_unloading_sim)**2))/len(porosity_unloading))
            k1_save=k1
            kp_save=kp/k1
            RMSE_data=[i,k1_save,RMSE_loading,kp_save,RMSE_unloading]
            RMSE_save.append(RMSE_data)
            if i<1:
                k1=k1+5000
                kp=kp
                continue
            if RMSE_loading>4000:
                k1=k1-(5000*((RMSE_save[i-1][3]-RMSE_save[i-2][3])/(RMSE_save[i-1][2]-RMSE_save[i-2][2])))
            else:
                k1=k1-(3000*((RMSE_save[i-1][3]-RMSE_save[i-2][3])/(RMSE_save[i-1][2]-RMSE_save[i-2][2])))
            kp=10*k1
            kc=0.1*k1
            ks=0.1*k1

    def checkdata():
        comp_pressure_save=abs(O.forces.f(plate.id)[2])/cross_area
        start_d=aabbExtrema()[0]+(0.2*D,0.2*D,0.2*D)
        end_d=aabbExtrema()[1]-(0.2*D,0.2*D,0.2*D)
        porosity_voxel=voxelPorosity(resolution=600,start=start_d,end=end_d)
        comp_data=[porosity_voxel,comp_pressure_save]
        compression_profile_save.append(comp_data)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hello,

> The code is not going through the loop

the "code" DOES go through the loop (as you describe yourself below) and does exactly what it is tell to do inside. Which is: set materials, particles, engines etc. and then directly go to the next for-loop iteration.

> The for loop runs by it self initially and even before simulation starts, and i set to 99 when the simulation starts.
> How can I run my simulation in the loop, for n number of iteration. And changing the parameter for each iteration.

so start the simulation before "you set to 99". Just use some running command inside the for loop:
###
for i in range (0,n_iter,1):
    # everything what you have now
    O.run(); O.wait() # or other stop condition [1]
    adjustParametersAccordingToResult()
###
Python reads and executes the code line-by-line.
Here, without any running, it just goes to the next iteration, as it is said.

Cheers
Jan

[1] https://yade-dem.org/doc/user.html#stop-conditions

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#2

Hi Jan,

It does not update parameter for in loop, the code thinks it is just one run. I am getting error message such as 'continue' not properly in loop. The simulation does not go through the loop and update parameter.

Best,
Mithu

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#3

Hi Jan,

Just realised, the word iteration could be misunderstood. By iteration I did not mean number iteration during simulation. I meant number of times, I will update my parameter.

Best,
Mithu

Revision history for this message
Jan Stránský (honzik) said :
#4

Well, of course, if you have "SyntaxError: 'continue' not properly in loop", it does nothing.
Python do not execute anything in case of syntax error.

First, try your simulation in one run, without any loop. You would have the same problem.
You have "def stopUnloading()" (starting a new scope), inside if condition and there "continue", which cannot be there (as the error says).

After you tune your single simulation, you can put the same code inside a for loop.

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#5

Hi Jan,

The main problem is that "code" DOES go through the loop, even if I remove the continue and update my parameter.

Best,
Mithu

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#6

The loop runs by itself first. Then the simulation starts.

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#7

Hi Jan,

I tried with simpler code, where I just update my parameter for each run. Where increase k1 with same value for each run.
It still doesn't go through the loop. i goes right upto 79.

Here is the code:

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d

o = Omega()

# Physical parameters
fr = 0.54
rho = 1550
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 50000
kp = 2.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.005
Cyl_height=0.045
cross_area=math.pi*(Tab_rad**2)

n_iter=80
RMSE_save=[]
#*************************************
for i in range (0,n_iter,1):
    o.reset()
    # Add material
    mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

    # Spheres for compression

    sp=pack.SpherePack()
    sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-18*Diameter),(3.0*Diameter,3.0*Diameter,18.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=527)
    sp.toSimulation(material=mat1)
    walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

    # Add engines
    o.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                            Bo1_Wall_Aabb(),
                            Bo1_Facet_Aabb()
                            ]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
        Ig2_Facet_Sphere_ScGeom(),
        Ig2_Wall_Sphere_ScGeom()],
        [Ip2_LudingMat_LudingMat_LudingPhys()],
        [Law2_ScGeom_LudingPhys_Basic()]
    ),
    NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
    PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
    #DeformControl(label="DefControl")
    ]

    def checkForce():
        if O.iter < 10000:
            return
        highSphere = 0.0
        for b in O.bodies:
            if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
                highSphere = b.state.pos[2]
            else:
                pass

        O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
        global plate
        plate = O.bodies[-1]
        plate.state.vel = (0, 0, -10)
        #O.engines = O.engines + [PyRunner(command='checkdata()', iterPeriod=1500)]
        fCheck.command = 'unloadPlate()'

    def unloadPlate():
        if abs(O.forces.f(plate.id)[2]) > Comp_force:
            plate.state.vel *= -1
            fCheck.command = 'stopUnloading()'

    def stopUnloading():
        if abs(O.forces.f(plate.id)[2]) == 0:
            o.wait()
            global k1
            global kp
            RMSE_data=[i,k1,kp]
            RMSE_save.append(RMSE_data)
            k1=k1+5000
            kp=i*k1
            kc=0.1*k1
            ks=0.1*k1

Revision history for this message
Jan Stránský (honzik) said :
#8

Hello,

all the problems are already explained in #1, please read it carefully.

in the OP:
> The code is not going through the loop
#5:
> The main problem is that "code" DOES go through the loop
#7:
> It still doesn't go through the loop.

Please be consistent.

> The loop runs by itself first. Then the simulation starts.

see #1, if you want to start the simulation before the entire loop is over, than just start simulation before the entire loop is over.

> I tried with simpler code, where I just update my parameter for each run. Where increase k1 with same value for each run.

No, you do not update parameters for each run, because there is NO RUN (see below or see #1). The posted code is actually expected to update nothing.

> It still doesn't go through the loop. i goes right upto 79.

what does "right upto" mean? If you have n_iter=80, it is logical that it loops up to 79 (range 80 is from 0 to 79).
79 does not appear from nowhere, in fact it is the result of perfect going through the (entire) loop.

My impression: the for loop is perfectly OK, doing exactly what it is told to do. I.e. set O.engines, add materials, add particles. Then, without any running, it jumps to the next iteration of the loop.
(Of course) it DOES NOT do what it is not told to do, i.e. it DOES NOT RUN anything (just because you have no start of the simulation in the for loop). See #1 how to solve it.

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#9

Hi Jan,

Sorry I was bit confused. I makes more sense now.
Correct if I am wrong, the code should be like this

Define initial values for the parameter(k1=1e6, kp=12*k1....)
for i in range (0,n_iter,1):
     o.reset()
     add spheres/material
     o.engine()
     o.run()
     update parameter based results()

Revision history for this message
Jan Stránský (honzik) said :
#10

Yes :-)

> O.run()

(not sure if it was just pseudo-code or meant to be used)
O.run() "as is" is not the right command, Python continues execution of next lines after O.run(), e.g. going to the next loop iteration and resetting the running simulation.

Either use
###
O.run(numIter,True) # True ... it waits until numIter is finished
###
or
###
O.run()
O.wait() # waits to O.pause() e.g. called from PyRunner checker, does not continue next line execution before simulation is paused
###

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#11

Hi Jan,

Thank you very much.
I can not open it, I am getting message: "script.py" is not responding and need Force Quit.
I am not sure if the code is wrong, or is due the capacity of my computer. Could you please have quick look and see if the code looks right.

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d

o = Omega()

# Physical parameters
fr = 0.54
rho = 1550
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 100000
kp = 2.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.005
Cyl_height=0.045
cross_area=math.pi*(Tab_rad**2)

Comp_press= 1e8
Comp_force=Comp_press*cross_area

n_iter=10
RMSE_save=[]
#*************************************
for i in range (3,n_iter,1):
    o.reset()
    # Add material
    mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

    # Spheres for compression

    sp=pack.SpherePack()
    sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-18*Diameter),(3.0*Diameter,3.0*Diameter,18.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=527)
    sp.toSimulation(material=mat1)
    walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

    # Add engines
    o.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
                            Bo1_Wall_Aabb(),
                            Bo1_Facet_Aabb()
                            ]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
        Ig2_Facet_Sphere_ScGeom(),
        Ig2_Wall_Sphere_ScGeom()],
        [Ip2_LudingMat_LudingMat_LudingPhys()],
        [Law2_ScGeom_LudingPhys_Basic()]
    ),
    NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
    PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
    #DeformControl(label="DefControl")
    ]

    def checkForce():
        if O.iter < 10000:
            return
        highSphere = 0.0
        for b in O.bodies:
            if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
                highSphere = b.state.pos[2]
            else:
                pass

        O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
        global plate
        plate = O.bodies[-1]
        plate.state.vel = (0, 0, -10)
        #O.engines = O.engines + [PyRunner(command='checkdata()', iterPeriod=1500)]
        fCheck.command = 'unloadPlate()'

    def unloadPlate():
        if abs(O.forces.f(plate.id)[2]) > Comp_force:
            plate.state.vel *= -1
            fCheck.command = 'stopUnloading()'

    def stopUnloading():
        if abs(O.forces.f(plate.id)[2]) == 0:
            o.pause()
            global k1
            global kp
            RMSE_data=[i,k1,kp]
            RMSE_save.append(RMSE_data)
            k1=k1+5000
            kp=i*k1
            kc=0.1*k1
            ks=0.1*k1
    O.run()
    O.wait()

Revision history for this message
Jan Stránský (honzik) said :
#12

> I can not open it

open what?

> I am getting message: "script.py" is not responding and need Force Quit.

where / how you get the message?

> I am not sure if the code is wrong, or is due the capacity of my computer. Could you please have quick look and see if the code looks right.

the code look s OK.
Try:
- run just one simulation without any loop, maybe it just takes very long
- run yade with "-n" option (no GUI, which should not be needed for this case)

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#13

>where / how you get the message?
When I run the the code in terminal, I am getting error message saying "myscript.py is not responding".

>run just one simulation without any loop, maybe it just takes very long
-The simulation does not work with one simulation.
-However if I remove o.run; o.wait from the code above. The simulation work for one simulation.

>run yade with "-n" option (no GUI, which should not be needed for this case)
-How do I run the with "-n" options

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#14

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

Revision history for this message
Jan Stránský (honzik) said :
#15

> When I run the the code in terminal, I am getting error message saying "myscript.py is not responding".

where you get the error? Int the terminal? From GUI?

> -The simulation does not work with one simulation.

Then it does not make sense to run it in a loop :-)

What does "does not work" mean?

> -However if I remove o.run; o.wait from the code above. The simulation work for one simulation.

What does "work" mean? Like what you do to make it work, e.g. click "play" button in GUI

> -How do I run the with "-n" options

yade -n script.py

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#16

Hi Jan,

Thanks, the codes runs. The issue is that it was extremely slow. Took upto 1h for one simulation.

Best,
Mithu

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#17

Hi Jan,

I am getting error from GUI.

Best,
Mithu

Revision history for this message
Jan Stránský (honzik) said :
#18

> I am getting error from GUI.

GUI is not responding while "O.wait"ing, so it is expected behavior.
For such "batch" simulations (either true yade-batch or this "O.reset in a for loop" approach), GUI should not be used.

Cheers
Jan

Revision history for this message
Mithushan Soundaranathan (mithushan93) said :
#19

Hi Jan,

Thank you very much :-)

Best,
Mithu