Saving/Loading

Hello,

> O.save(namefile) saves engines, interactions, ..., but it does not save parameters or other things.
> But it seems that all the history is simply disappearing, for example, volumetric strains.

What are "...". "parameters", "other things", "the history", "volumetric strains" etc.?

Please provide a MWE [1] showing what you mean.

> I am obliged to save (O.save) close Yade and then later on open it and load what I have saved.

Just for our information, why and what is the use case?

> My question is: Is there any technique to save the history of volumetric and axial strains?

Definitely yes, there are many options actually.
Which to use depends on what you want to achieve.

> Hope my message was clear.

not really, sorry..

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Hey Jan,
Sorry for not being very clear.

 > What are "...". "parameters", "other things", "the history", "volumetric strains" etc.?

Basically, I am trying to couple DEM with LBM in order to model partially saturated medias. And LBM is written in C++, and it's the main code in which I call Yade. For example, for some number of LBM iterations (n_LBM), I launch Yade for (n_DEM) iterations. At the end of the Yade cycle (at n_DEM iteration), I save two things:
- I save some physical quantities in a text file, such as volumetric strain, stresses, axial strains, kinetic energy, etc. In each row, I have the number of iterations and the physical parameters mentioned above.
- I save the simulation by using O.save(filename). (Indeed, for loading the simulation again in Yade, I use O.load(filename))
Once Yade is finished and everything is saved, Yade exists and the main code (LBM) has the control again, and so on. This the DEM-LBM coupling algorithm.

The example I am currently studying is Triaxial test in partially saturated media. Indeed, in order to verify that the algorithm works well, I have simulated a dry REV, which means that I commented the LBM block code so that only Yade is used but using my algorithm (cycles of opening and closing Yade). And I compare the results of my algorithm's simulation with a normal Yade simulation.
In both simulations, sigma_zz is increasing (in the saved text file), which means that the upper wall is compressing my REV while maintaining sigma_xx = sigma_yy = sigma_conf (which is a good sign). However, using my algorithm, epsilon_zz, epsilon_vv, and kinetic energy are all zeros while using a normal Yade triaxial simulation, epsilon_zz, epsilon_vv, and kinetic energy are changing.
Somehow, when I load the simulation using O.load() strains and kinetic energy are reset.
Hope my message is more precise now.

Best
Nabil

Thanks for some clarification

> Once Yade is finished ..., Yade exists

Yade exists? Yade exits?

> Hope my message is more precise now.

More, but not fully..
Please provide a MWE [1]

Cheers
Jan

Hi (again),

TriaxialStressController.strain is not O.saved if that is your concern.

Nevertheless, it will be re-computed at most TriaxialStressController.computeStressStrainInterval iterations after the O.load, I think, to a correct value since TriaxialStressController.*0 (e.g. TriaxialStressController.depth0) are themselves saved.

As such, I'm not expecting real issues for the strain quantities but you can tell us the contrary with a MWE ;-)

For the kinetic energy, it would also help if you tell us exactly how you're measuring it (which YADE / Python commands).

(I would bet bodies velocities are correctly O.saved and O.loaded)

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

Hello again,

I don't know how to add an attachment, but I will try to explain the obtained results using the coupling that I have already written in the previous answer. All I want to verify is that the coupling is working. Therefore, I will suppose that all capillary forces will be equal to 0 (as if I am studying a dry sample) and I will compare this simulation with a normal YADE simulation without calling it from a C++ code.

I have tested the coupling for 300000 iterations for YADE, it means that the C++ code will open YADE and run it for 300000 iterations. After these iterations, YADE will exit, C++ will take control again and then launch YADE for another 300000 iterations.

From 0 to 300000 iterations, the volumetric strain is increasing (in absolute value, I am making a triaxial test). The volumetric strain at iteration number 300000 is -0.000124, the sigma_zz = -103356.446315, and void ratio is = 0.747134.

Next, YADE exits, and then will be opened via C++, and then O.load() is applied to load the particles and etc.. Therefore, YADE will start from 300001 iteration and finishes at iteration number 600001 to exit and reopen again.

It is interesting to note that at iteration number 300001 (just after re-opening YADE), the volumetric strain becomes -2E-06 (approx 0), sigma_zz =-103362.362661, and void ratio is 0.747133.

So there is a continuity for stresses, but it seems not working for the strains, and this is exactly what I was asking for. How to make the strain continuous when reloading data?

Hope that this time is more clear, thank you so much in advance!
Nabil
..

Hi (again) Nabil,

Trying to rephrasing myself, it is impossible (unless you modify C++ source code, obviously) to get back correct TriaxialStressController.strain just after a O.load.

Playing with TriaxialStressController.computeStressStrainInterval (possibly while the model is imposed to freeze, see b.state.blockedDOFs, if you wish to) may diminish the inconvenience, which is just about post-processing in my view.