Removing particles during thermal cycling

Asked by Zoheir Khademian on 2020-12-14

Hello All,

I would like to remove a particle from an assembly after a few Hydro-thermo-mechanical cycles (using thermal and flow engines) and then continue cycling the model. After removing the particle, I enforce re-triangulation [1] and then use [2] but I get this error and Yade crushes once I continue running the model:

python3.6: /usr/include/boost/smart_ptr/shared_ptr.hpp:734: typename boost::detail::sp_member_access<T>::type boost::shared_ptr<T>::operator->() const [with T = yade::Body; typename boost::detail::sp_member_access<T>::type = yade::Body*]: Assertion `px != 0' failed.
Aborted (core dumped)

I am wondering what I am missing here and would appreciate your help.

Thanks

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngineT.updateTriangulation
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.emulateAction

Here is the MWE:

from yade import pack, ymport, plot, utils, export, timing
import numpy as np

young=5e6
mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)
O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600e10,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600e10,label='spheres'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=0.333,num=200,seed=11)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 NewtonIntegrator(damping=0.5)
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1

tri_pressure = 1000
triax.goal1=triax.goal2=triax.goal3=-tri_pressure
triax.stressMask=7
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.1 and abs(-tri_pressure-triax.meanStress)/tri_pressure<0.001:
    break

triax.internalCompaction=False

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 10
flow.meshUpdateInterval=1
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,10,10]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]
flow.tZero=25

flow.dead=0
thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=True
thermal.solidThermoMech = True
thermal.fluidThermoMech = True
thermal.advection=True
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,0,1]
thermal.thermalBndCondValue=[0,0,0,0,0,45]
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2700
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0

timing.reset()
O.dt=1e-4
O.dynDt=False
thermal.dead=0
flow.emulateAction()

def bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos(0.025,0.025,0.025)

from yade import plot

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature((0.025,0.025,0.025)),
  p=flow.getPorePressure((0.025,0.025,0.025)),
  t=O.time,
  i = O.iter,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))}

plot.plot(subPlots=False)
O.run(5,1)

#Removing one particle
Removedbody = bodyByPos(0.,0.,0.)
O.bodies.erase(Removedbody.id)
updateTriangulation=1
O.run()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Zoheir Khademian
Solved:
2021-01-18
Last query:
2021-01-18
Last reply:
2020-12-14
Robert Caulk (rcaulk) said : #1

Hey Zoheir,

Thanks for reporting this with a working MWE. I think this is now fixed, if you want to test before it is merged you can clone and compile the branch titled "eraseBodyThermalEngine" [1], otherwise you will need to wait for it to work its way into yadedaily which could take a couple weeks.

Cheers,

Robert

[1]https://gitlab.com/yade-dev/trunk/-/tree/eraseBodyThermalEngine

Robert Caulk (rcaulk) said : #2

Also a side note: your line:

updateTriangulation = 1

Does nothing inside yade. I will let you brainstorm why that is the case and how to fix it.

Zoheir Khademian (zkhademian) said : #3

Robert, thank you so much for the quick reply. I truly appreciate it. I will recompile it and report back to you after trying the MWE.

> updateTriangulation = 1 Does nothing inside yade. I will let you brainstorm why that is the case and how to fix it.

Is it because I did not associate it with the engine? like: "flow.updateTriangulation=1"?
If so, Is "flow.updateTriangulation=1" the same as [1], using flow.meshUpdateInterval=1 after removing the particle and then return it to flow.meshUpdateInterval=-1?

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.meshUpdateInterval

Robert Caulk (rcaulk) said : #4

>Is it because I did not associate it with the engine? like: "flow.updateTriangulation=1"?

Yes :-)

>If so, Is "flow.updateTriangulation=1" the same as [1], using flow.meshUpdateInterval=1 after removing the particle and then return it to flow.meshUpdateInterval=-1?

flow.updateTriangulation is a boolean, so it is a True or False (but 1 or 0 work too). If it is set to True, the next iteration will force flowEngine to retriangulate. Then it will set updateTriangulation back to false.

meshUpdateInterval stays set at the number you choose, and will force flowEngine to update the triangulation every *meshUpdateInterval* iterations. This will need to stay at 1 if you are using the full THM coupling.

Cheers,

Robert

Robert Caulk (rcaulk) said : #5

This should be in an answered state, user must have accidentally clicked "still need help."

Zoheir Khademian (zkhademian) said : #6

This solved my problem