yade-users team mailing list archive

[Anqi H <question683336@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Anqi H is still having a problem:
Hi Deepak,

Thank you so much for your help so far. I've modified my script based on
your comment to check if a body does not have any cohesive bonds in its
interactions, if true then delete this sphere from fluidcoupling IDs. In
the first openFoam time step, some spheres were deleted as expected,
however, after the deletions and in the same foam timestep, it seems the
program hangs after 500 yade iterations and in the fluidCoupling step (I
have put the particle deletion pyrunner before calling fluidCoupling). I
have copied my python script below. I am sorry it's probably difficult
to run this script as it requires two other txt files that stored the
sphere locations. I was wondering if it was due to some obvious errors
in my script that I haven't noticed?


from __future__ import print_function
import sys
from yadeimport import *
from yade.utils import *
from yade import ymport

initMPI()                           #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling();     #Initialize the engine
fluidCoupling.getRank();            #part of Initialization.


#example of spheres in shear flow : two-way point force coupling
class simulation():

  def __init__(self):

    O.periodic = True;
    
    #proppant properties
    FrictAng_p = 0.9
    Density_p = 2650
    Young_p = 100e6
    TensileStr_p=3000
    Cohesion_p=3000

    Young = 57e8 #apparent modulus
    FrictAng = 0.5
    Density = 2650
    Poisson = 0.28
    Cohesion = 38e6 # pa
    TensileStr = 38e6 # pa

    rock = JCFpmMat(type=1,young=Young,frictionAngle=FrictAng,density=Density,poisson=Poisson,tensileStrength=TensileStr,cohesion=Cohesion,label='rock')
    proppant = JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')
    O.materials.append(JCFpmMat(young=Young_p,frictionAngle=0,density=0,label='wallmat'))

    O.materials.append(proppant)
    O.materials.append(rock)

    rock_assembly =
O.bodies.append(ymport.textExt('new_rock.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=rock))

    for b in O.bodies:
      if b.state.pos[1]>0.017 or b.state.pos[1]<0.0135:
        O.bodies.erase(b.id)
      #if (b.state.pos[1]>0.0135 and b.state.pos[1]<0.015) :
       # b.groupMask=2
      else:
        b.dynamic = False
        b.state.vel=(0, 0, 0)

    proppant_assembly = O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
    print ('length of bodies proppant '+str(len(O.bodies)))

    for b in proppant_assembly:
      O.bodies[b].groupMask=2
      if O.bodies[b].state.pos[0] < 0:
        print("found it ")
        O.bodies.erase(b)


    sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere and b.material.label=='proppant')]

    #coupling engine settings
    fluidCoupling.setNumParticles(len(sphereIDs))
    fluidCoupling.setIdList(sphereIDs)
    fluidCoupling.isGaussianInterp=False;  #use pimpleFoamYade for gaussianInterp

    # Integrator
    newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
     # add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.


    O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2, label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2, label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
  ),
  GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
        PyRunner(command='sim.deletePar()',iterPeriod=30, label='checkPar'),
        fluidCoupling, #to be called after timestepper
        PyRunner(command='sim.printMessage()', iterPeriod= 1000, label='outputMessage'),
  newton,
        #PyRunner(command='sim.deletePar()',iterPeriod=50, label='checkPar'),
        VTKRecorder(fileName='yadep/3d-vtk-',mask = 2,recorders=['spheres','colors'],iterPeriod=1000)
    ]
      

  def printMessage(self):

     print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
     if O.iter == 4000:
         maxVel = 0.05
         for b in O.bodies:
             if type(b.shape)==Sphere:
                 bodyVel = abs(b.state.vel.norm())
                 if bodyVel > maxVel:
                     raise ValueError("Body velocity exceeds imposed shear velocity by ", abs(bodyVel-maxVel))

  def deletePar(self):
    ids = fluidCoupling.getIdList()
    for b in ids:
      temp=0
      for i in O.bodies[b].intrs():
        if i.phys.isCohesive==True:
          temp+=1
      if temp==0:
        fluidCoupling.eraseId(b)
        print('deleted '+str(b))
    print("******YADE-ITER = " + str(O.iter) +" **********")

    #if O.bodies[b].state.nbInitBonds==0 or
(O.bodies[b].state.nbInitBonds==O.bodies[b].state.nbBrokenBonds):


  def steprun(self):
    O.step()
    bols.aabbEnlargeFactor = 1
    ig2s.interactionDetectionFactor = 1
    print('-------------------reset aabbEnlargeFactor------------------')

  def irun(self,num):
      O.run(num,1)


if __name__=="__main__":
    sim = simulation()
    sim.steprun()
    sim.irun(5000) #run 5000 iteration and wait
    # print("body id = ", O.bodies[34].id)
    fluidCoupling.killMPI()

import builtins
builtins.sim=sim

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