New question #700611 on Yade: https://answers.launchpad.net/yade/+question/700611
Hi all, I am running the script [1] in OpenMpi with increased domain and smaller particles. [1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/FluidCouplingPFV/drainage-2PFV-Yuan_and_Chareyre_2017.py The pack is filled with 680 particles of radius 5e-4, 1520 particles of radius 3e-4 and 3590 particles of radius 2e-4. This pack (matrix_vtest.txt) was inported to the drainage simulation 2PFV [1] and each particle was substituted by agglomarate of particles as: particle size 5e-4 were substituted by agglomerates of particles of radius5e-5 (agg5e4_5e5.txt) particle size 3e-4 were substituted by agglomerates of particles of radius 3e-5 (agg3e4_3e5.txt) particle size 2e-4 were substituted by agglomerates of particles of radius 2e-5 (agg2e4_2e5.txt) When I run the drainage simulation 2PFV [1] without substituting the pack by agglomerate (just to check its functionality) it works. But, when I substitute the particles by the agglomerates I get the following error when the drainage starts: double free or corruption (!prev) [node030:278860] *** Process received signal *** [node030:278860] Signal: Aborted (6) [node030:278860] Signal code: (-6) [node030:278860] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x138e0)[0x2ae56996a8e0] [node030:278860] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x141)[0x2ae5699b6e71] [node030:278860] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x112)[0x2ae5699a0536] [node030:278860] [ 3] /lib/x86_64-linux-gnu/libc.so.6(+0x7e2b8)[0x2ae5699f82b8] [node030:278860] [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x85d0a)[0x2ae5699ffd0a] [node030:278860] [ 5] /lib/x86_64-linux-gnu/libc.so.6(+0x8757c)[0x2ae569a0157c] [node030:278860] [ 6] /usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN4yade3CGT7NetworkINS0_12_TesselationINS0_18TriangulationTypesINS_18TwoPhaseVertexInfoENS_16TwoPhaseCellInfoEEEEEE19defineFictiousCellsEv+0x1c5)[0x2ae56d60ae55] [node030:278860] [ 7] /usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN4yade38TemplateFlowEngine_TwoPhaseFlowEngineTINS_16TwoPhaseCellInfoENS_18TwoPhaseVertexInfoENS_3CGT12_TesselationINS3_18TriangulationTypesIS2_S1_EEEENS3_25FlowBoundingSphereLinSolvIS7_NS3_18FlowBoundingSphereIS7_EEEEE18buildTriangulationEdRSB_+0xa18)[0x2ae56d6341b8] [node030:278860] [ 8] /usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN4yade18TwoPhaseFlowEngine14initializationEv+0x4f)[0x2ae56d5a4e8f] [node030:278860] [ 9] /usr/lib/x86_64-linux-gnu/yade/libyade.so(_ZN5boost6python7objects23caller_py_function_implINS0_6detail6callerIMN4yade18TwoPhaseFlowEngineEFvvENS0_21default_call_policiesENS_3mpl7vector2IvRS6_EEEEEclEP7_objectSH_+0x43)[0x2ae56d5aff03] [node030:278860] [10] /usr/lib/x86_64-linux-gnu/libboost_python39.so.1.74.0(_ZNK5boost6python7objects8function4callEP7_objectS4_+0x2bb)[0x2ae56e8596db] [node030:278860] [11] /usr/lib/x86_64-linux-gnu/libboost_python39.so.1.74.0(+0x20968)[0x2ae56e859968] [node030:278860] [12] /usr/lib/x86_64-linux-gnu/libboost_python39.so.1.74.0(_ZN5boost6python21handle_exception_implENS_9function0IvEE+0x73)[0x2ae56e85e893] [node030:278860] [13] /usr/lib/x86_64-linux-gnu/libboost_python39.so.1.74.0(+0x1e0c2)[0x2ae56e8570c2] [node030:278860] [14] /usr/bin/python3.9(_PyObject_MakeTpCall+0x39b)[0x51df8b] [node030:278860] [15] /usr/bin/python3.9[0x53c575] [node030:278860] [16] /usr/bin/python3.9(_PyEval_EvalFrameDefault+0x5418)[0x517508] [node030:278860] [17] /usr/bin/python3.9[0x510d5d] [node030:278860] [18] /usr/bin/python3.9(_PyEval_EvalCodeWithName+0x47)[0x510b07] [node030:278860] [19] /usr/bin/python3.9(PyEval_EvalCode+0x23)[0x5f38b3] [node030:278860] [20] /usr/bin/python3.9[0x5f8540] [node030:278860] [21] /usr/bin/python3.9[0x529ee4] [node030:278860] [22] /usr/bin/python3.9(_PyEval_EvalFrameDefault+0x52b)[0x51261b] [node030:278860] [23] /usr/bin/python3.9[0x510d5d] [node030:278860] [24] /usr/bin/python3.9(_PyFunction_Vectorcall+0x342)[0x529362] [node030:278860] [25] /usr/bin/python3.9(_PyEval_EvalFrameDefault+0x52b)[0x51261b] [node030:278860] [26] /usr/bin/python3.9[0x511657] [node030:278860] [27] /usr/bin/python3.9(_PyFunction_Vectorcall+0x342)[0x529362] [node030:278860] [28] /usr/bin/python3.9(_PyEval_EvalFrameDefault+0x52b)[0x51261b] [node030:278860] [29] /usr/bin/python3.9[0x511657] [node030:278860] *** End of error message *** I wounder if it has to do with Numpy/Scipy not supporting such a large number of in the sparce matrix and if would have a work around. Or if it is something else. Below is the script I modified to import the pack and agglomerates (.txt). I did not provide an external link to these txt files, but I can do that if necessary. #!/usr/bin/python # -*- encoding=utf-8 -*- import os from yade import mpy as mp from yade import pack from yade import bodiesHandling from yade import export from yade import utils from yade import ymport #import math ############################################ ### DEFINING VARIABLES AND MATERIALS ### ############################################ # The following 5 lines will be used later for batch execution nRead=readParamsFromTable( num_spheres=3000,# number of spheres compFricDegree = 1, # contact friction during the confining phase key='_triax_base_', # put you simulation's name here unknownOk=True ) from yade.params import table num_spheres=table.num_spheres# number of spheres key=table.key targetPorosity = 0.40 #the porosity we want for the packing compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process) finalFricDegree = 30 # contact friction during the deviatoric loading rate=0 # loading rate (strain rate) damp=0.8 # damping coefficient stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below) #2e4+70e4medio 1e4+70e4bom 1e4+60e4bom 3e4+90e4+w3,1,-1-the best young=20e5 # contact stiffness200e4 young2=20e5 youngcoat=20e5 bondstr=0.3e3 bondstr2=0.3e3 bondstrcoat=1e3 ## create materials for spheres and plates mat=O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres')) O.materials.append(JCFpmMat(type=1,young=20e7,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls')) O.materials.append(JCFpmMat(type=1,young=youngcoat,poisson=0.3,frictionAngle=radians(1),density=1500,tensileStrength=bondstrcoat,cohesion=bondstrcoat,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstrcoat,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spherescoat')) ## create walls around the packing mn,mx=Vector3(0,0,0),Vector3(0.01005,0.01005,0.01005) mnbox,mxbox=Vector3(-0.0001,-0.0001,-0.0001),Vector3(0.01005,0.0115,0.01005) walls=aabbWalls([mnbox,mxbox],thickness=0,material='walls') wallIds=O.bodies.append(walls) O.bodies.append(ymport.textExt("matrix_vtest.txt", format='x_y_z_r', shift=Vector3(0,0,0), scale=1.0,material='spherescoat',color=(0,1,1))) ################Particle substitution by large aggregate###################################################################### bodid=[] a=[] for b in O.bodies: if b and isinstance(b.shape,Sphere): # print (b.shape.radius) if b.shape.radius==0.0005: bodid.append(b.id) a.append(b.state.pos) i=0 for p in bodid: t=a[i] f1=O.bodies.append(ymport.textExt("agg5e4_5e5.txt", format='x_y_z_r', shift=t-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1))) O.bodies.erase(bodid[i]) i=i+1 bodidd=[] aa=[] for bb in O.bodies:# in sp: if bb and isinstance(bb.shape,Sphere): # print (bb.shape.radius) if bb.shape.radius==0.0003: bodidd.append(bb.id) aa.append(bb.state.pos) ii=0 for pp in bodidd: tt=aa[ii] f2=O.bodies.append(ymport.textExt("agg3e4_3e5.txt", format='x_y_z_r', shift=tt-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1))) O.bodies.erase(bodidd[ii]) ii=ii+1 bodiddd=[] aaa=[] for bbb in O.bodies:# in sp: if bbb and isinstance(bbb.shape,Sphere): print (bbb.shape.radius) if bbb.shape.radius==0.0002: bodiddd.append(bbb.id) aaa.append(bbb.state.pos) iii=0 for ppp in bodiddd: ttt=aaa[iii] f3=O.bodies.append(ymport.textExt("agg2e4_2e5.txt", format='x_y_z_r', shift=ttt-Vector3(0,0,0), scale=1.0,material='spheres',color=(0,1,1))) O.bodies.erase(bodiddd[iii]) iii=iii+1 ############################################################################################################################## #====================================================================================================== triax=TriaxialStressController( stressMask = 0, internalCompaction=False, # If true the confining pressure is generated by growing particles wall_front_activated=True, wall_back_activated=True, wall_top_activated=True, wall_bottom_activated=True, wall_left_activated=True, wall_right_activated=True, goal1=-5, goal2=-15, goal3=-5, ) newton=NewtonIntegrator(damping=damp) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1,label='Physicspheres')],#,xSectionWeibullShapeParameter=1.5, xSectionWeibullScaleParameter=1 [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3), triax, newton ] ################################################### ### REACHING A SPECIFIED POROSITY PRECISELY ### ################################################### while triax.porosity>targetPorosity: setContactFriction(radians(compFricDegree)) print ("\r Friction: ",compFricDegree," porosity:",triax.porosity) # porosity will decrease O.run(100,1) print ("### Compacted state saved ###") print(triax.stress(3)[1]) # Change contact friction (remember that decreasing it would generate instantaneous instabilities) setContactFriction(radians(finalFricDegree)) #======================================= #set stress control on x and z, we will impose strain rate on y triax.stressMask = 2 triax.wall_bottom_activated=0 #now goal2 is the target strain rate triax.goal1=rate#triax.stress(1)[0] triax.goal3=rate#triax.stress(5)[2] triax.goal2=triax.stress(3)[1] ##################################################### ### Example of how to record and plot data ### ##################################################### #from yade import plot from yade import plot O.run(100,True) #strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L (strain) ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2] si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2] ## a function saving variables def history(): plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2, ev=-triax.strain[0]-triax.strain[1]-triax.strain[2], s11=-triax.stress(triax.wall_right_id)[0]-si0, s22=-triax.stress(triax.wall_top_id)[1]-si1, s33=-triax.stress(triax.wall_front_id)[2]-si2, e=math.exp(-triax.strain[1]-ei1)-1, pc=-unsat.bndCondValue[2], sw=unsat.getSaturation(False), z1=O.bodies[3].state.pos[1], i=O.iter) O.engines=O.engines+[PyRunner(iterPeriod=1000,command='history()',label='recorder')] ####################################################### ## Drainage Test under oedometer conditions ### ####################################################### ##Instantiate a two-phase engine unsat=TwoPhaseFlowEngine() ##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant unsat.bndCondIsPressure=[0,0,1,1,0,0] unsat.bndCondIsWaterReservoir=[0,0,1,0,0,0] unsat.bndCondValue=[0,0,-1e8,0,0,0] unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir unsat.initialization() unsat.surfaceTension = 0.0728 #start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt f5=open('SwPc5.txt',"w") ts=O.dt pgstep= 4500000*ts #30Pa/s print (pgstep) pgmax= 3000#9316 #Pa for pg in arange(1.0e-8,pgmax,pgstep): unsat.bndCondValue=[0,0,(-1.0)*pg,0,0,0] unsat.invasion() unsat.computeCapillaryForce() unsat.meshUpdateInterval=500 unsat.defTolerance=-1 unsat.updateTriangulation=True print(unsat.getSaturation(False),pg,-triax.strain[1]) for b in O.bodies: O.forces.setPermF(b.id, unsat.fluidForce(b.id)) while 1: O.run(100,True) unb=unbalancedForce() if unb<0.05: break f5.write(str(pg)+" "+str(unsat.getSaturation(False))+" "+str(triax.strain[1])+"\n") f5.close() mp.DOMAIN_DECOMPOSITION= True #automatic splitting/domain decomposition mp.mpirun(1) #passive mode run -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
