yade-users team mailing list archive

[Mithushan Soundaranathan <question697379@xxxxxxxxxxxxxxxxxxxxx>]

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

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)





        


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