THM slow for large number of subdomains

[jvwilliams23]

Hi all,

I am running an old THM script very similar to those in #25555. What I am finding, is that for a large number of channels/junctions (even when number of elements per channel is very low), the solve is very slow (initialisation is also a bit slow, but it is not prohibitive). Do you think there is any potential fix for this? @GiudGiud @joshuahansel

See the two test cases:

1. 100 channels, each with n_elem=32 (3299 elements and 199 subdomains). Setting up time = 4.9 s. Execution time = 5.16 s.
   straight-pipe-100.i.txt

2. 800 channels each with n_elem=4 (3999 elements and 1599 subdomains). Setting up time = 318.65 s. Execution time = 296.3 s. straight-pipe-800.i.txt

3. (extra test) 100 channels with n_elems=100 (100099 elements and 199 subdomains). Setting up time = 21 s. Execution time = 53 s.

I am profiling it now, I will post the output below.

Thanks,
Josh

test-2 log:

*** Info ***

The following occurred in the AuxKernel 'pipe-0-1:area_aux' of type ProjectionAux.

Projection lowers order, please expect a loss of accuracy
Still Setting Up............................................................             [318.65 s] [ 4613 MB]

Framework Information:
MOOSE Version:           git commit 78413b32bf on 2026-01-29
LibMesh Version:         502d42871efb6093844d2e02abe82bb15d48a033
PETSc Version:           3.24.3
SLEPc Version:           3.24.0
Current Time:            Thu Feb  5 14:24:43 2026
Executable Timestamp:    Tue Feb  3 13:35:28 2026

Input File(s):
  /work4/scd/scarf1324/project/farscapeWS4/chimera/chimera-coolant-newprof/test/straight-pipe-800.i

Command Line Argument(s):
  -w

Command Line Input Argument(s):
  GlobalParams/n_elems=4

Checkpoint:
  Wall Time Interval:      Every 3600 s
  User Checkpoint:         Disabled
  # Checkpoints Kept:      2
  Execute On:              TIMESTEP_END

Parallelism:
  Num Processors:          1
  Num Threads:             1

Mesh:
  Parallel Type:           replicated
  Mesh Dimension:          3
  Spatial Dimension:       3
  Nodes:                   4799
  Elems:                   3999
  Num Subdomains:          1599

Nonlinear System:
  Num DOFs:                13595
  Num Local DOFs:          13595
  Variables:               { "rhoA" "rhoEA" "rhouA" } { "rhoV" "rhouV" "rhovV" "rhowV" "rhoEV" }
  Finite Element Types:    "MONOMIAL" "MONOMIAL"
  Approximation Orders:    "CONSTANT" "CONSTANT"

Auxiliary System:
  Num DOFs:                41597
  Num Local DOFs:          41597
  Variables:               "A" "A_linear" { "H" "P_hf" } "T" "e" "p" { "rho" "v" } "vel" { "vel_x" "vel_y"
                             "vel_z" }
  Finite Element Types:    "MONOMIAL" "LAGRANGE" "MONOMIAL" "MONOMIAL" "MONOMIAL" "MONOMIAL" "MONOMIAL"
                             "MONOMIAL" "MONOMIAL"
  Approximation Orders:    "CONSTANT" "FIRST" "CONSTANT" "CONSTANT" "CONSTANT" "CONSTANT" "CONSTANT"
                             "CONSTANT" "CONSTANT"

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


Currently Executing
  Performing Initial Setup
    Setting Up Initial Conditions.                                                       [ 10.91 s] [ 5545 MB]
    Setting Up Materials
      Computing Initial Material Values                                                  [ 12.91 s] [ 5545 MB]
    Finished Setting Up Materials                                                        [ 12.95 s] [ 5545 MB]
  Still Performing Initial Setup.
    Computing User Objects.........

                                                      [ 52.98 s] [ 5377 MB]
  Finished Performing Initial Setup                                                      [ 96.35 s] [ 5377 MB]


Time Step 0, time = 0

Time Step 1, time = 0.0001, dt = 0.0001


    Computing Residual
      Computing User Objects........                                                     [ 46.56 s] [ 5380 MB]
    Still Computing Residual                                                             [ 57.73 s] [ 5380 MB]
    Computing User Objects.....                                                          [ 36.64 s] [ 5380 MB]

 0 Nonlinear |R| = 1.046590e+01

    Computing Jacobian
      Computing User Objects.......                                                      [ 46.21 s] [ 5380 MB]
    Still Computing Jacobian                                                             [ 52.55 s] [ 5380 MB]
    Computing Residual                                                                   [ 13.44 s] [ 5380 MB]

 1 Nonlinear |R| = 9.355692e-05

    Computing Jacobian                                                                   [ 10.66 s] [ 5380 MB]
    Computing Residual                                                                   [ 13.51 s] [ 5380 MB]

 2 Nonlinear |R| = 1.088543e-06

 Solve Converged!

  Finished Solving                                                                       [192.26 s] [ 5380 MB]


Outlier Variable Residual Norms:
  rhoEA: 7.420131e-07
  rhoEV: 7.964555e-07

Steady-State Relative Differential Norm: 8.89291e-05

Finished Executing                                                                       [296.30 s] [ 5380 MB]

[jvwilliams23]

Profiler output

pprof_out.txt

[joshuahansel]

Hi fellow Josh!

I am running an old THM script

So are you using an old executable too? We made changes to the volume junction since then so that it doesn't introduce a new scalar variable for each junction.

[jvwilliams23]

Hi Josh 😄! No I compiled latest MOOSE a couple of days ago.

[joshuahansel]

I did two runs, one with your 800-channel case as-is, and one where I replaced the volume junctions with one-to-one junctions:

• Volume junctions: nodes=4799,elems=3999,subdomains=1599,dofs=13595: 186.6 s total
• 1-to-1 junctions: nodes=4000,elems=3200,subdomains=800,dofs=9600: 89.8 s total

The volume junction changes probably solved one scaling issue and created another: now volume junctions share common variable names, but each create a new block (NodeElem), so they double the number of subdomains here.

Thank you for the profile - that's very helpful! @GiudGiud So if I'm reading this right, we're spending a fair amount of time checking if the subdomain has changed? I'm trying to determine how much we could possibly trim off here. 127 s (AttribSubdomainsMatch) of the 1146 s (total)?

[jvwilliams23]

Thanks for checking, just an FYI this is a simpler version of a case where we need to use VolumeJunctions as the junctions have more than one connection.

[joshuahansel]

Understood, I was trying to see if I could definitively determine that it was the number of subdomains, but of course then I realized those volume junctions have extra DoFs on them, so it's not a perfect comparison. Regardless, I think that you've demonstrated a clear trend. I think it's just a question of if there is anything reasonable to do about it.