# yade-users team mailing list archive

## [Question #697379]: Optimise material parameter

```New question #697379 on Yade:

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
import pandas as pd
import numpy as np
from PIL import Image
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)

Cyl_height=0.045

Comp_force=Comp_press*cross_area

n_iter=80
RMSE_save=[]
#*************************************
for i in range (0,n_iter,1):
o.reset()
compression_profile_save=[]
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)

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)]

if abs(O.forces.f(plate.id)[2]) > Comp_force:
plate.state.vel *= -1

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)
k1_save=k1
kp_save=kp/k1
RMSE_save.append(RMSE_data)
if i<1:
k1=k1+5000
kp=kp
continue
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)

--