Problems about the flow.multithread

I've been using the "flow.multithread" to accelerate the calculation speed of problems about saturated soil recently. I've encounter several problems about the "flow.multithread".

1. When I write "multithread=True" before the simulation begins, just like the following, it reports that Segmentation fault (core dumped).

########
O.engines=[
 ...
 FlowEngine(multithread=1, dead=0, label="flow", ompThreads=10),
 ...
]
########

However, if I don't set multithread=True at the first time I define the flow engine, set flow.multithread=True after running several steps instead, it works and the calculation is accelerated. I have no idea why this happens.

2. When I set flow.multithread=True after running several steps, although it works, there appears to be a period error. After certain steps, all pore pressure turn to zero and the next step everything is back to normal again. This error occurs at a fixed interval.

I'm using the
Linux version 5.4.0-47-generic (buildd@lgw01-amd64-038) (gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04))
Yade 20201007-4243~57d1c9c~bionic1

Here's a MWE to produce the error. Through the VTK file can get the period error.

######################################
from yade import pack,export

SoilMat=CohFrictMat(young=20e9,poisson=0.3,density=2650,frictionAngle=0.7,alphaKr=50,alphaKtw=50,momentRotationLaw=True,label='soil')
O.materials.append((SoilMat))

mn,mx=Vector3(-.25,-.25,0),Vector3(.25,.25,.5)
walls=aabbWalls([mn,mx],thickness=0,material='soil')
wallIds=O.bodies.append(walls)

snum=1000
sp=pack.SpherePack()
sp.makeCloud((-.25,-.25,0),(.25,.25,.5),rMean=.001,rRelFuzz=.3,num=snum)
sp.toSimulation(material=SoilMat)

Target_Stress=-1000000

triax=TriaxialStressController(
 stressMask = 7,
 thickness = 0,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

Pressure=10000
flow=FlowEngine(
 dead=1,
 defTolerance=0.001,
 meshUpdateInterval=200,
 useSolver=3,
 #multithread=True,
 permeabilityFactor=0.001,
 bndCondIsPressure=[0,0,0,0,1,1],
 bndCondValue=[0,0,0,0,0,Pressure], #The depth of the water is 5.1m
 boundaryUseMaxMin=[0,0,0,0,0,0],
 viscosity=0.00089 #The dynamic viscosity of water is 8.90*10-4 Pa`s at about 25 degree.
)

newton=NewtonIntegrator(damping=0.1)

#Engine of the simulation
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 flow,
 triax,
 PyRunner(command="print('%s: Unbalncedforce %s, Meanstress %s' % (O.iter,unbalancedForce(),triax.meanStress))",realPeriod=1,label='Monitor'),
 newton,
 PyRunner(command="flow.saveVtk()",iterPeriod=100,dead=1,label='vtksaver2'),
]

triax.goal1=triax.goal2=triax.goal3=Target_Stress

while 1:
  O.run(1000,True)
  if unbalancedForce()<0.1 and abs((Target_Stress-triax.meanStress)/Target_Stress)<0.01:
    break

Monitor.dead=1
vtksaver2.dead=0
flow.dead=0
triax.internalCompaction=False
triax.stressMask=0
triax.goal1=triax.goal2=triax.goal3=0
O.run(1,1)
flow.multithread=True
######################################