yade-users team mailing list archive

[Nabil YOUNES <question702176@xxxxxxxxxxxxxxxxxxxxx>]

Question #702176 on Yade changed:
https://answers.launchpad.net/yade/+question/702176

Nabil YOUNES posted a new comment:
Hello everyone,
Sorry for replying a bit late..
Here's my MWE:

C/C++ the LBM code in which capillary forces are computed, then YADE is called for n_DEM iterations: 
// Code C/C++

Block code 1

// This will call the triax_test.py script (shown below)

system((std::string("yade -n -x triax_test.py ") + std::to_string(npDEM)
+ std::string(" ") +  std::to_string(currentLoad_DEM)).c_str());

Block code 2
When the code reads the command system((std::string("yade ....)), Yade will open and it will take control:

YADE:

# ==================================== #
from yade import plot
import os
import sys
import time
import numpy as np

##############################################################################
# Specific functions
##############################################################################

def check():
    """
    Check the convergence of each loading step of the triaxial test and swich to the next
    step when converged.
    """
    if O.iter % 10000 == 0 and O.iter != 0:
        global currentLoad, isoStrain, loadingName    
        print('I AM IN CHECK AT ITERATION ', O.iter)
        print('\n')
        unb=unbalancedForce()
        if currentLoad==0:
            mean=triax.meanStress
    #        test=unb<0.001 and \
            test=unb<0.0001 and \
                abs(mean-loadPath[currentLoad][0])/loadPath[currentLoad][0]<0.01
            print('At iter %i, mean=%0.3e, unb=%0.3e and loadPath=%0.3e'%(O.iter,mean,unb,loadPath[currentLoad][0]))
            if test:
                isoStrain[0]=triax.strain[0]
                isoStrain[1]=triax.strain[1]
                isoStrain[2]=triax.strain[2]
                currentLoad=1
                triax.stressMask=mask[currentLoad]
                triax.goal1=loadPath[currentLoad][0]
                triax.goal2=loadPath[currentLoad][1]
                triax.goal3=loadPath[currentLoad][2]
                triax.max_vel=10 # In order to allow the lateral wall to move fast enough
                # recorder.dead=False
                # dataSaver.dead=False
                # sampleSaver.dead=False
                
        elif currentLoad==1:
            test=abs(triax.strain[2]-isoStrain[2])>0.15
            q=triax.stress(triax.wall_back_id)[2]-triax.stress(triax.wall_left_id)[0]
            p=(triax.stress(triax.wall_left_id)[0]+\
            triax.stress(triax.wall_bottom_id)[1]+\
            triax.stress(triax.wall_back_id)[2])/3
            eta=q/p
            print('At iter %i, axial stain reached %0.3e and eta=%0.3f'%(O.iter,triax.strain[2]-isoStrain[2],eta))
            if test:
                O.save('endTriax_'+loadingName+'.yade.gz')
                plot.saveDataTxt('triax_'+loadingName+'_data.txt')
            #    O.pause()
        X_particles_max = max([b.state.pos[0] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
        Z_particles_max = max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
        Y_particles_max = max([b.state.pos[1] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
        X_particles_min = max([b.state.pos[0] - b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
        Z_particles_min = max([b.state.pos[2] - b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])       
        Y_particles_min = max([b.state.pos[1] - b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])       
        if X_particles_max > 3e-3 or Y_particles_max > 3e-3 or Z_particles_max > 3e-3 or X_particles_min < -3e-3 or Y_particles_min < -3e-3 or Z_particles_min < -3e-3:
            print("URGENT !!!! Boders have been crossed\n\n")
            O.pause()
            O.save('endTriax_'+loadingName+'.yade.gz')
            plot.saveDataTxt('triax_'+loadingName+'_data.txt')

        eps_zz=triax.strain[2]-isoStrain[2]
        if abs(eps_zz) > 0.3:
            print("PAUSE BEACUSE epz_zz > 0.3 \n\n")
            O.pause()


# save_data_triax() function that saves multiple parameters
def save_data_triax():
    if O.iter % 500 == 0 and O.iter != 0:
        if not os.path.isfile('triax_data.dat'):
            f = open("triax_data.dat", "w")
            f.write('iter time[sec] KineticEnergy ElastPotential PlastDissip sig_xx sig_yy sig_zz Eps_xx Eps_yy Eps_zz VoidIndex Eta Eps_vv\n')
            f.close()
        else:
            f = open("triax_data.dat", "a")
            q=triax.stress(triax.wall_back_id)[2]-triax.stress(triax.wall_left_id)[0]
            p=(triax.stress(triax.wall_left_id)[0] + triax.stress(triax.wall_bottom_id)[1] + triax.stress(triax.wall_back_id)[2])/3
            eta=q/p
            # Stresses
            sig_xx=triax.stress(0)[0]
            sig_yy=triax.stress(2)[1]
            sig_zz=triax.stress(4)[2]
            # Strains
            eps_xx=triax.strain[0]-isoStrain[0]
            eps_yy=triax.strain[1]-isoStrain[1]
            eps_zz=triax.strain[2]-isoStrain[2]
            # Void Index
            e=porosity()/(1-porosity())
            f.write('%d %f %f %f %f %f %f %f %f %f %f %f %f %f \n' % (O.iter, O.iter * O.dt, O.energy['kinetic'], O.energy['elastPotential'], O.energy['plastDissip'], sig_xx, sig_yy, sig_zz, eps_xx, eps_yy, eps_zz, e, eta, eps_xx + eps_yy + eps_zz))
            f.close()


def Save_data():
	print("\n SAVING DATA \n")
#	t0_Saving_Data = time.clock() 
#	O.save('data.xml.bz2')
#	tf_Saving_Data = time.clock() - t0_Saving_Data
#	print("Time of saving data.xml.bz2 is ", tf_Saving_Data)
	t0_Saving_Data = time.clock() 
	O.save('data.gz')
	tf_Saving_Data = time.clock() - t0_Saving_Data
	print("Time of saving data.gz is ", tf_Saving_Data)
	f = open("LB_data.dat", "w+")
	f.write('id r(m) x(m) y(m) z(m) vx(m/s) vy(m/s) vz(m/s)\n')
	f.write("HAVE_YADE_PAUSED?: %f \n" % test)
	f.write("currentLoad %f \n" % currentLoad)
	for b in O.bodies:
		if type(b.shape)==Sphere:
			f.write("%d %0.10f %0.10f %0.10f %0.10f %0.10f %0.10f %0.10f \n" % (b.id, b.shape.radius,b.state.pos[0],b.state.pos[1], b.state.pos[2],b.state.vel[0],b.state.vel[1], b.state.vel[2]))
	f.close()


O.dt = 1e-7


loadingName='looseSamplebis'
filename = 'TextureTest'
sigmaIso=-1*10**5
mask=[7,3]
strainSpeed=0.01
loadPath=[[sigmaIso,sigmaIso,sigmaIso],[sigmaIso,sigmaIso,-strainSpeed]]
nbLoads=len(mask)
#currentLoad has the values of 0 for confining the sample and 1 for axial strain while sigma0 is fixed
currentLoad = int(sys.argv[2])
test = 0


# To store the isotropic stain level
isoStrain=[0,0,0]

# Deviatoric states to save
#strainTarget=[0.01,0.02,0.03,0.04]
etaTarget=[0.,0.2,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65]
EpsTarget=[0.1]
progress=0
# to keep in memory the number of saved texture
nbTextureSaved=0

# reload the sample from .py file
young=100e6
poisson=0.42
frictAngle=35
density=3000
damp=0.05

O.load('data.gz')

# Enable energy tracking
O.trackEnergy=True


# Stiffness of walls and spheressotropic initial compression
triax=TriaxialStressController(label='triax',
                # Isotropic initial compression
                stressMask=mask[currentLoad],
                goal1=loadPath[currentLoad][0],
                goal2=loadPath[currentLoad][1],
                goal3=loadPath[currentLoad][2],
                internalCompaction=False,                
                # Comment if internalCompaction=True
                max_vel=0.01,stressDamping=0.25,strainDamping=0.8,
                # Uncomment if internalCompaction=True
                maxMultiplier=1.0001,
                finalMaxMultiplier=1.000001
                )
                
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
        InteractionLoop(
                 [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
                 [Ip2_FrictMat_FrictMat_FrictPhys()],
                 [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
#        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,
#                                   timestepSafetyCoefficient=0.8),
        triax,      
        NewtonIntegrator(damping=damp,label='newton')
        ]

checker=PyRunner(iterPeriod=1,command='check()',label='checker')

O.engines+=[checker]

O.run(int(sys.argv[1]) + 1, True)

Save_data()

# ==================================== #

Once YADE exits, the command line (in C++) will go to Block code 2.

Hope it's better now. Thank you so much in advance.
Nabil

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