temp=Nan?

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

My purpose is to achieve high-temperature heating of the rock mass to simulate a high-temperature environment, so I wrote the following script.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.