About the examples/ThermalEngine

Hello,

>>1.When I run the original script (unchanged), the problem is the same as before: plot is blank and body.State.Temp = Nan.

Will you please install the latest source code version? You appear to have installed a snapshot from mid-March.

>>So I guess the problem has something to do with the number of spheres?

Sphere size and distribution changes the maximum allowable timestep for numerical stability. If you decrease sphere size, you will need to decrease O.dt to maintain stability.

Cheers,

Robert

Increasing sphere number in the same geometrical domain size, means smaller spheres. Smaller spheres means smaller characteristic lengths for the system. Smaller char. lengths means smaller allowable timestep. Please refer to this discussion for more details:

https://answers.launchpad.net/yade/+question/698948

Cheers,

Robert

Thanks a lot Robert!

>>Will you please install the latest source code version? You appear to have installed a snapshot from mid-March.

Yes,I also noticed that.As you suggested[1], I used the following command:

cd trunk
git checkout -b working_state 61c90697
cd ../build
cmake -DCMAKE_INSTALL_PREFIX=../install ../trunk
make -j4 install

And then the version changes to 2021-03-16.
I compiled it from the source code again. The version is yade-2021-09-28.git-52f0c5a.But when I run the original script (unchanged),the same problem still exists...

>>Sphere size and distribution changes the maximum allowable timestep for numerical stability. If you decrease sphere size, you will need to decrease O.dt to maintain stability.

Ok,I will decrease O.dt and provide timely feedback.

Hello Robert!
I have the latest findings, which may help solve this problem:

I follow your suggestion and consider changing O.dt and the number of spheres respectively.I found something interesting..

1.I have tried the number of spheres as 100,200,300,500,1000,2000,and the O.dt as 1e-4 to 1e-6 (I also tried using the O.dt=PWaveTimeStep(), not sure if it's appropriate).The success of the simulation does not seem to have much to do with O.dt but O.dt affects the speed of the simulation(The premise is that the simulation is successful, that is, there is temperature instead of Nan)
2.As described in 1, I found that when the number of spheres is 100,200,1000,2000, the simulation is successful (i.e. the pore and body have their own temperature rather than Nan).But 300,500,the problem still exists,like the original script (unchanged).With these scripts,there is a same warning:
CHOLMOD warning: matrix not positive definite. file: ../Supernodal/t_cholmod_super_numeric.c line: 911
(Maybe this is the reason of problem?)
3.As mentioned in 2, even in the successful cases, the curves are different.In the example script, I guess the ideal curve should be to keep the pore temperature at 343.15 and gradually increase the particle temperature from 333.15 to 343.15.If I'm right, this is true only when num = 100, and in other cases, the body temperature will increase to 353.15..
What is wrong with the problem?

Best regards.

Hello,

>>I compiled it from the source code again. The version is yade-2021-09-28.git-52f0c5a.But when I run the original script (unchanged),the same problem still exists...

I cannot reproduce this issue using latest source. The example scripts run as intended. It may be an issue with distribution, you are running 20.04. Can you please install latest yade source on Ubuntu 18.04 and confirm your issue?

>>I have tried the number of spheres as 100,200,300,500,1000,2000,and the O.dt as 1e-4 to 1e-6 (I also tried using the O.dt=PWaveTimeStep(), not sure if it's appropriate).The success of the simulation does not seem to have much to do with O.dt but O.dt affects the speed of the simulation(The premise is that the simulation is successful, that is, there is temperature instead of Nan)

Have you actually looked at your packing? It seems you are simulating a gaseous state of spheres (although there is no way to know without an MWE). ThermalEngine was not designed to handle this kind of packing.

Please provide an MWE.

Cheers,

Robert

Hello,Robert

1.I have followed your suggestion, installed ubuntu18.04 and installed Source-code.Unfortunately,the problem still exists..:
################
CHOLMOD warning: matrix not positive definite. file: ../Supernodal/t_cholmod_super_numeric.c line: 911
/usr/lib/python3/dist-packages/matplotlib/__init__.py:831: MatplotlibDeprecationWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  mplDeprecation)
/usr/lib/python3/dist-packages/matplotlib/__init__.py:801: MatplotlibDeprecationWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  mplDeprecation)
################

printAllVersions():
Yade version : 2021-10-08.git-a6d7547
Yade features : LOGGER USEFUL_ERRORS VTK OPENMP GTS QT5 CGAL PFVFLOW PFVFLOW LINSOLV MPI TWOPHASEFLOW FEMLIKE GL2PS LBMFLOW THERMAL PARTIALSAT POTENTIAL_PARTICLES POTENTIAL_BLOCKS
Yade config dir: ~/.yade-2021-10-08.git-a6d7547
Yade precision : 53 bits, 15 decimal places, without mpmath, PrecisionDouble
```

Libraries used :

| library | cmake | C++ |
| ------------- | -------------------- | ------------------- |
| boost | 106501 | 1.65.1 |
| cgal | | 4.11 |
| clp | 1.16.11 | 1.16.11 |
| cmake | 3.10.2 | |
| coinutils | 2.10.14 | 2.10.14 |
| compiler | /usr/bin/c++ 7.5.0 | gcc 7.5.0 |
| eigen | 3.3.4 | 3.3.4 |
| freeglut | 2.8.1 | |
| gl | | 20190911 |
| ipython | 5.5.0 | |
| metis | | 5.1.0 |
| mpi | 3.1 | ompi:2.1.1 |
| mpi4py | 2.0.0 | |
| mpmath | 1.0.0 | |
| openblas | | OpenBLAS 0.2.20 |
| python | 3.6.9 | 3.6.9 |
| qglviewer | | 2.6.3 |
| qt | | 5.9.5 |
| sphinx | 1.6.7-final-0 | |
| sqlite | | 3.22.0 |
| suitesparse | 5.1.2 | 5.1.2 |
| vtk | 6.3.0 | 6.3.0 |

```
Linux version : Ubuntu 18.04.6 LTS
Architecture : amd64
Little endian : True

(ps:other example scripts such as flowscenario and THMcoupling can run normally,only this script is unsuccessful, maybe it is my personal installation problem?)

2.>>Have you actually looked at your packing? It seems you are simulating a gaseous state of spheres (although there is no way to know without an MWE). ThermalEngine was not designed to handle this kind of packing.

Sorry for my stupidity,does the gaseous state of the sphere mean that the sphere is too scattered (not dense enough?).I refer to thermoHydroMechanical_coupling.py before using makecloud. Is there anything wrong?

MWE as follows:
#################################
from yade import pack, ymport
from yade import timing
import numpy as np
import shutil
timeStr = time.strftime('%m-%d-%Y')
num_spheres=1000# number of spheres
young=1e9
rad=0.003

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05) # corners of the initial packing

thermalCond = 2. #W/(mK)
heatCap = 710. #J(kg K)
t0 = 333.15 #K

r = rad
k = 2.0 # 2*k*r
Cp = 710.
rho = 2600.
D = 2.*r
m = 4./3.*np.pi*r**2/rho

identifier = '-noFlowScenario'

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)

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(3),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=SpherePack()
sp.makeCloud(mn,mx,num=100,rRelFuzz=0.333)
sp.toSimulation(material='spheres')

print('num bodies ', len(O.bodies))

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

ThermalEngine = ThermalEngine(dead=1,label='thermal');

newton=NewtonIntegrator(damping=0.2)
intRadius = 1
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow",multithread=False),#introduced as a dead engine for the moment, see 2nd section
 ThermalEngine, 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
]

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

flow.dead=0
flow.defTolerance=-1
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=343.15
flow.pZero=0

thermal.dead=0
thermal.debug=False
thermal.fluidConduction=True
thermal.ignoreFictiousConduction=True
thermal.conduction=True
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
thermal.bndCondIsTemperature=[0,0,0,0,0,0]
thermal.thermalBndCondValue=[0,0,0,0,0,0]
thermal.fluidK = 0.6069
thermal.fluidConductionAreaFactor=1.
thermal.uniformReynolds=10 # using an extremely low reynolds number to maintain a heat transfer coefficient
thermal.particleT0 = 333.15
thermal.particleDensity=2600.
thermal.particleK = 2. #thermalCond
thermal.particleCp = 710.
thermal.tsSafetyFactor=0
thermal.useKernMethod=True
thermal.useHertzMethod=False
timing.reset()

O.dt=0.1e-4
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)

#print "found body of interest at", bodyOfInterest.state.pos

axis = np.linspace(mn[0], mx[0], num=5)
axisBodies = [None] * len(axis)
axisTrue = np.zeros(len(axis))
for i,x in enumerate(axis):
 axisBodies[i] = bodyByPos(x, mx[1]/2, mx[2]/2)
 axisTrue[i] = axisBodies[i].state.pos[0]

from yade import plot

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature((0.025,0.025,0.025)),
  #ftemp2=flow.getPoreTemperature((0.5,0.1,0.5)),
  #ftemp3=flow.getPoreTemperature((0.5,0.9,0.5)),
  p=flow.getPorePressure((0.025,0.025,0.025)),
  t=O.time,
  i = O.iter,
  temp1 = axisBodies[0].state.temp,
  temp2 = axisBodies[1].state.temp,
  temp3 = axisBodies[2].state.temp,
  temp4 = axisBodies[3].state.temp,
  temp5 = axisBodies[4].state.temp,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp

)
 plot.saveDataTxt('txt'+timeStr+identifier+'/temps'+identifier+'.txt',vars=('t','i','p','ftemp1', 'temp1','temp2','temp3','bodyOfIntTemp'))

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

def endFlux():
 if O.time >= 30:
  O.pause()
O.engines=O.engines+[PyRunner(iterPeriod=10,command='endFlux()')]
from yade import plot

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

#print "starting thermal sim"
O.run()
###########################
Thanks Robert again!

Hello,

Please confirm that you have latest source compiled on 18.04 and noFlowScenario.py (no modifications) does not run.

If you confirm this, then there is something else going on which may be difficult to trouble shoot. I am running 18.04 on latest source and noFlowScenario.py runs as expected.

>>does the gaseous state of the sphere mean that the sphere is too scattered

Yes it means you have no spheres touching one another, which would imply quite large conduction diffusivity, which would imply very small timestep. But besides this, the conduction model is not physically realistic in such a gaseous state so it should not be used.

You can view your packing by clicking "Show 3D" in the GUI window.

Cheers,

Robert

Hello,
I am sorry to tell you that I do use latest source compiled on 18.04(This can be verified by the information in printallversions() ).
The only change is that:
shutil.copyfile(sys.argv[0],'txt'+timeStr+identifier+'/'+sys.argv[0])
(FileNotFoundError: [Errno 2] No such file or directory,So I delete this line,but this should not lead to problems?)
In fact,I guess the root of the problem lies in 'CHOLMOD warning'?(Just guess)

>>Yes it means you have no spheres touching one another, ....
I rechecked thermoHydroMechanical_coupling.py and found that after makecloud, triaxial internal compaction is carried out to make the sphere in close contact.Thank you for your explanation. I have realized my mistake.

(ps:I seem to remember that during the compilation and installation process (i.e. during the execution of the 'make' command), the terminal prompts warnings about Law2_ScGeom_MindlinPhys_Mindlin?,does it have anything to do with this)

Best regards.

Hello,

noFlowScenario.py works as expected on 18.04 using latest source code.

I ran the script you provided, it also works.

Since I am unable to reproduce the error that you describe, I cannot help you.

Cheers,

Robert

Hello,

It seems that the situation is very bad for me personally...
Ok,I will try to solve this problem myself..

Finally, how should CHOLMOD warning be solved.

Thank you for your timely and patient help again!

Just an update:
Hello everyone,I seem to have solved this problem by changing the value of flow. Usesolver(Although I'm not sure if it's the right thing to do)

Best regards.

Hello,

Although you have not specified what value you changed flow.useSolver to, I can tell you that anything besides 4 will not solve for advection. Hence why it is important to make sure the example scripts run as is.

Cheers,

Robert

Hello Robert,
The situation is like this, I consulted the forum and referred to the solution of [1] to change the value of flow.useSolver.
As you said,when I change the value to 0 or 3,it seems that the error not happen and everything is normal.
I also refer to the explanation of flow.useSolver in [2]
(ps:Maybe it’s because I don’t know what the exact result should look like)

Yes,the script have thermal.advection=True for achieving advection,but when I run the script as useSolver=4,it still have the same problem although I have update the yadedaily to 20211018-5981~9b2e38e.:
CHOLMOD warning: matrix not positive definite. file: ../Supernodal/t_cholmod_super_numeric.c line: 911
the plot is blank and O.bodies.state.temp=nan..
What should I do.

Best regards.

[1]https://answers.launchpad.net/yade/+question/697642
[2]https://yade-dem.org/doc/yade.wrapper.html?highlight=flowengine#yade.wrapper.FlowEngine.useSolver

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Hello,

Please confirm that updating your yade installation has solved your problem. If not, please indicate the exact changes you have made compared to the example script, and which change in particular is the one that results in your issue.

Thank you,

Robert

Avoiding launchpad janitor

This question was expired because it remained in the 'Needs information' state without activity for the last 15 days.

Marking as answered to avoid launchpad janitor until OP returns.

That solved my problem.

Hello,

Please either indicate which answer solved your problem by marking it with a star, or describe the solution to your problem.

Thank you for helping keep an organized knowledge base.

-rc