temp=Nan using FlowEngine and ThermalEngine

Dear all,

I use yadedaily:
Yade 20210909-5909~5b11526~focal1
Using python version: 3.8.10 (default, Jun 2 2021, 10:49:15)
[GCC 9.4.0]
Ubuntu 20.04

Since my ultimate goal is to achieve high-temperature rock testing, I decided to apply high-temperature conditions to the sample first.However,in the process of my attempts, there are several questions related to FlowEngine and ThermalEngine that confused me.
1.Since the purpose of my script is to achieve high temperature conditions, I don’t think I need to use FlowEngine.But when I run the script,an error message appears: 'the core has been dumped'(Ps:Chinese on my system).When I added FlowEngine, this error disappeared. Does it mean that ThermalEngine cannot be used alone?
2.As mentioned above, if two engines must be used together, then what is the difference between FlowEngine.thermalBndCondValue and ThermalEngine.thermalBndCondValue?
3.The script I ended up running is as follows.Unfortunately, it was not as successful as I expected.I found that the body.state.temp=Nan.What's the problem?

##########code############
from yade import pack,ymport,plot
from yade import timing
import numpy as np
import shutil
timeStr=time.strftime('%m-%d-%Y')
mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)
identifier = '-keep_temp'
young=1e9

if not os.path.exists('VTK'+timeStr+identifier):
 os.mkdir('VTK'+timeStr+identifier)
else:
 shutil.rmtree('VTK'+timeStr+identifier)
 os.mkdir('VTK'+timeStr+identifier)
if not os.path.exists('txt'+timeStr+identifier):
 os.mkdir('txt'+timeStr+identifier)
else:
 shutil.rmtree('txt'+timeStr+identifier)
 os.mkdir('txt'+timeStr+identifier)

O.materials.append(FrictMat(young=1e9,poisson=0.25,frictionAngle=radians(3),density=2600,label='spheres'))
O.materials.append(FrictMat(young=1e9,poisson=0.25,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=O.bodies.append(ymport.text('radius_0.0008.sphere',color=(0.1,0.1,0.9),material='spheres'))
print('num bodies:',len(O.bodies))

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

Thermal = ThermalEngine(dead=1,label='thermal');
newton=NewtonIntegrator(damping=0.2)
intRadius=1.5
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow",multithread=False),
 Thermal,
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 VTKRecorder(iterPeriod=500,fileName='VTK'+timeStr+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=1,label='VTKrec'),
 newton
]

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

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.dynamic=False

flow.dead=0
flow.defTolerance=-1 #0.3
flow.meshUpdateInterval=-1
flow.useSolver=4
flow.permeabilityFactor= 1
flow.viscosity= 0.001
flow.bndCondIsPressure=[0,0,0,0,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.debug=False
flow.fluidRho = 997
flow.fluidCp = 4181.7
flow.bndCondIsTemperature=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]
flow.tZero=298.15
flow.pZero=0

thermal.dead=0
thermal.debug=False
thermal.ignoreFictiousConduction=True
thermal.conduction=True
thermal.bndCondIsTemperature=[1,1,1,1,1,1]
thermal.thermalBndCondValue=[473.15,473.15,473.15,473.15,473.15,473.15]
thermal.particleT0=298.15
thermal.particleDensity=2600
thermal.particleK=2
thermal.particleCp=710
thermal.useKernMethod=True
thermal.uniformReynolds=10
thermal.useKernMethod=True

timing.reset()

O.dt=0.1e-3
O.dynDt=False

O.run(1,1)
flow.dead=0

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
 print('found closest body ', cBody.id, ' at ', cBody.state.pos)
 return cBody

bodyOfInterest = bodyByPos(0.025,0.025,0.025)
TolTemp=523.15

def Stop_condition():
 if bodyOfInterest.state.temp > TolTemp :
  O.stop()

O.engines=O.engines+[PyRunner(iterPeriod=100,command='Stop_condition()')]

def history():
 plot.addData(
  t=O.time,
  i=O.iter,
  bodyOfIntTemp=O.bodies[bodyOfInterest.id].state.temp
  )
 plot.saveDataTxt('txt'+timeStr+identifier+'/temps'+identifier+'.txt',vars=('t','i','bodyOfIntTemp'))

O.engines=O.engines+[PyRunner(iterPeriod=500,command='history()',label='recorder')]
VTKrec.dead=0

plot.plots={'i':('bodyOfIntTemp')}
plot.plot()
O.saveTmp()

print('starting thermal sim')
O.run()
#############################

Many thanks for your help!
Best regards.