Thermal hydraulic model long initialisation times

[jvwilliams23]

Hi,

I am running a thermal hydraulic analysis of a coolant system. The graph has about 808 FlowChannels and 557 VolumeJunctions to connect them.

The issue I am having is when I begin to run the simulation (using the 'combined' module solver) it takes a long time at the setup stage of "Initializing Equation Systems". To do this, it takes 5000 seconds and uses 814 MB memory. The mesh only has 3232 nodes and 2424 elements. I have noticed this time increases when I use more elements and also more VolumeJunctions (if I use JunctionOneToOne1Phase for all junctions with 1 inlet and 1 outlet, it takes like 300 seconds to initialise, which is still quite long).

Is there any way to avoid this? It seems there is a variable created for every "junction" (also these show up in paraview as individual variables), which must cause the issue. I have attached a log file below.

Best wishes,
Josh

--------------------------------------------
Setting Up
  Finished Setting Up Options                                                            [  5.94 s] [  186 MB]


*** Info ***
Projection lowers order, please expect a loss of accuracy
Still Setting Up......
  Initializing
    Initializing Equation Systems........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................... [5027.27 s] [  814 MB]
  Finished Initializing                                                                  [5030.18 s] [ 2804 MB]
Finished Setting Up                                                                      [5079.07 s] [ 2341 MB]

Framework Information:
MOOSE Version:           git commit a643d739e2 on 2023-06-26
LibMesh Version:
PETSc Version:           3.16.6
SLEPc Version:           3.16.2
Current Time:            Fri Sep 22 11:07:09 2023
Executable Timestamp:    Tue Sep 12 15:36:10 2023

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh:
  Parallel Type:           replicated
  Mesh Dimension:          1
  Spatial Dimension:       3
  Nodes:                   3232
  Elems:                   2424
  Num Subdomains:          808

Nonlinear System:
  Num DOFs:                10057
  Num Local DOFs:          10057
  Variables:               { "junction-0:rhoEV" "junction-0:rhoV" "junction-0:rhouV" "junction-0:rhovV"
                             "junction-0:rhowV" ... "junction-99:rhoEV" "junction-99:rhoV" "junction-99:rhouV"
                             "junction-99:rhovV" "junction-99:rhowV" } { "rhoA" "rhoEA" "rhouA" }
  Finite Element Types:    "SCALAR" "MONOMIAL"
  Approximation Orders:    "FIRST" "CONSTANT"

Auxiliary System:
  Num DOFs:                34799
  Num Local DOFs:          34799
  Variables:               "A" "A_linear" { "H" "P_hf" "T" } "T_wall" "e" { "junction-0:T" "junction-0:p"
                             "junction-0:vel" "junction-100:T" "junction-100:p" ... "junction-97:p" "junction-97:vel"
                             "junction-99:T" "junction-99:p" "junction-99:vel" } { "p" "rho" "v" "vel_x"
                             "vel_y" "vel_z" }
  Finite Element Types:    "MONOMIAL" "LAGRANGE" "MONOMIAL" "LAGRANGE" "MONOMIAL" "SCALAR" "MONOMIAL"

  Approximation Orders:    "CONSTANT" "FIRST" "CONSTANT" "FIRST" "CONSTANT" "FIRST" "CONSTANT"

Execution Information:
  Executioner:             Transient
  TimeStepper:             ConstantDT
  TimeIntegrator:          ImplicitEuler
  Solver Mode:             NEWTON
  MOOSE Preconditioner:    SMP (auto)

[GiudGiud]

Hello

@joshuahansel will have more insights but I think the prime way to accelerate these things is to move from having one object per flow channel and junction to a single vectorized object adding all these components. This will take a little bit of reworking the code, and we typically only do that if we have a need for a faster simulation.

There might also be trivial slow downs of the code. If you profile the code we may just be able to fix them. Please see this page for profiling techniques:
https://mooseframework.inl.gov/application_development/profiling.html

Guillaume

[jvwilliams23]

I can provide an input file in few hours.

[jvwilliams23]

Here are two input files. One has 100 edges, 99 junctions and takes approx 70 seconds to initialise on my laptop. The other has 150 edges and 149 junctions, takes approx 300 seconds to initialise on my laptop.

straight-pipe-100.i.txt
straight-pipe-150.i.txt

[roystgnr]

Do we have any input files that replicate the problem?

This was just poor wording on my part - it was silly for me to digress about wanting to replicate a different problem in this ticket, not after we determined that that problem wouldn't get in the way of fixing this one.

On the other hand, posting input files for this problem was also a good idea; that'll make it much easier to test potential fixes. And such atrocious scaling... N elements, each trying to futz with its own O(N^2) dense matrix of SCALARs ... With O(N^3) total we'd expect a 1.5x larger problem to be 3.375x longer instead of 1.5x longer to initialize here, which is bad enough (it would imply that 100 junctions is already initializing 10,000x slower than it should be), but your numbers look even worse than that! I hope there's not something going on here I'm still missing...

I'm tempted to tackle this on two fronts... getting rid of the SCALAR use here as per #25593 is definitely the right thing to do in THM, but maybe it wouldn't hurt to improve the SCALAR handling in libMesh too. My ideas for that would at least get the THM scaling down to O(N^2), which is still much worse than the proper O(N) but which might be a sufficient quick fix for these small cases. Any SCALAR fix would eventually be mooted by the fixes at the THM level, but there are likely other users with large numbers of SCALAR variables who really do need them, so optimization there wouldn't go to waste.

[joshuahansel]

Ok I've been working on this idea of using field variables on the junctions. I tried it out on your problem with 3 different methods:

• Current: what we're doing currently, i.e., using scalar variables for junctions.
• New A: using field variables on NodeElems, keeping the same variable names unique to each junction, e.g., junction1:pressure.
• New B: using field variables on NodeElems, using common variable names for all junctions, e.g., pressure, which is fine because these are block-restricted.

Here are rough timing results in seconds, with the 100-edge case first and 150-edge case in parentheses, and the "?" where I didn't bother to try:

• Current: 110(380) setup, 116(390) total
• New A: 29(?) setup, 240(?) total
• New B: ~0(~0) setup, 5.4(9.8) total

So basically, if you just change to field variables but don't take advantage of the fact that you don't need unique variable names anymore, then you do significantly decrease the setup time, but for some reason, you increase the solve time (including residual and Jacobian computation) significantly, giving a worse total time. However, if you do start using common variable names, then you see that setup time go to basically zero (I didn't really measure it), while your solve time stays very good.

So this is great news. I don't yet have an answer to why "New A" is so terrible. I will need to actually use a profiler to do this (and I do plan to do this - maybe it will shed light on other performance pitfalls). But whatever the reason, THM has a good solution here, and I'll work on getting this merged soon.

Thank you to @roystgnr and @lindsayad for the critical suggestions, and thank you @jvwilliams23 for raising the issue and providing good input and examples! I'm hopeful that THM's performance for large problems will be greatly improved! 📈 🎉

[jvwilliams23]

Sounds great @joshuahansel. Could you point me to your branch with this code? I would like to try it out

[joshuahansel]

Ok I started #28514 as a WIP. I'm just finishing up fixing tests.

Note that it might require some prerequisite work: the new variable names in the system can cause different ordering of the variables in the nonlinear system, which can lead to convergence issues (#28503). All in all, I'm hopeful that everything should be merged by the end of next week 🤞

[joshuahansel]

@licharlot You have some experience with large numbers of components. We've already applied that dependency resolution fix, right?

[jvwilliams23]

Thanks! For one-to-one junctions I will just use a higher friction factor on the adjacent pipes as a form loss for now (will move to FormLossFromFunction1Phase later).
So the number of junctions is now approximately 68, and initialisation takes 100 seconds.

When I increase the number of elements per flowchannel from 3 to 20, it takes around 600 seconds.

[joshuahansel]

And to be clear, JunctionParallelChannels1Phase has a form loss option, so that can be used for any case where the connected channels are parallel.

[NorthMagic]

I have a similar problem with THM initialization out of time especailly with multi-cpu cores/threads parallel computing. I tried to simulate a full core with hundreds of parallel 1D flow channels and without junctions, it takes more than 300 s to finish the initialization. Moreover, I found that the initialization time increased greatly with the number of used cpu cores and threads, which was unacceptable for doing simulations.
From my understanding to the mentioned answers above, a long-time initialization was cuased by too many variables which did not exist in my simulation. More impotantly, the initialization time increased with used cpu cores. Could you please give me some suggestions?

[joshuahansel]

Ok so you're saying your simulation does not have junctions, so you shouldn't have the same issue. Can you also profile your simulation? https://mooseframework.inl.gov/application_development/profiling.html

[NorthMagic]

Sorry, I have never used this profiling tool. I need more time to install and profile my simulation. I put some new testing outputs in Discussion #27941 and the input file for running THM was posted.