Fail to converge in an eigen problem

[js-jixu]

Hi, Moltres developers.

I want to solve the steady state neutron diffusion equation with Moltres to get parameters like neutron flux distribution and fission heat. My model is a single rod model with fuel, clading and coolant. The keff I calculated with OpenMC is probably between 1.3-1.4. I wrote a new input file referring to the input file of moltres/problems/2021-cnrs-benchmark/phase-0/nts.i:

[GlobalParams]
  num_groups = 6
  num_precursor_groups = 6
  use_exp_form = false
  group_fluxes = 'group1 group2 group3 group4 group5 group6'
  pre_concs = 'pre1 pre2 pre3 pre4 pre5 pre6'
  temperature = 600
  sss2_input = true
  account_delayed = true
[../]

[Mesh]
  file = 3layers_3d_4parts.msh
[]

[Problem]
  type = FEProblem
[]

[Nt]
  var_name_base = group
  vacuum_boundaries = 'cool_bottom cool_top clad_bottom clad_top fuel_bottom fuel_top'
  create_temperature_var = false
  transient = false
  eigen = true
  scaling = 1e3
  block = 'fuel cool clad'
[]

[Precursors]
  [./pres]
    var_name_base = pre
    outlet_boundaries = ''
    constant_velocity_values = true
    u_def = 0
    v_def = 0
    w_def = 0
    nt_exp_form = false
    family = MONOMIAL
    order = CONSTANT
    loop_precursors = false
    transient = false
    eigen = true
    scaling = 1e3
    block = 'fuel cool clad'
  [../]
[]

[Materials]
  [./fuel]
    type = GenericMoltresMaterial
    property_tables_root = './xs/fuel_'
    interp_type = 'linear'
    block = 'fuel'
  [../]
  [./clad]
    type = GenericMoltresMaterial
    property_tables_root = './xs/clad_'
    interp_type = 'linear'
    block = 'clad'
  [../]
  [./cool]
    type = GenericMoltresMaterial
    property_tables_root = './xs/cool_'
    interp_type = 'linear'
    block = 'cool'
  [../]
[]

[Executioner]
  type = InversePowerMethod
  max_power_iterations = 100

  # fission power normalization settings
  normalization = 'powernorm'
  normal_factor = 21230           # Watts, 280e6 / (109 * 11 * 11)

  xdiff = 'group1diff'
  bx_norm = 'bnorm'
  k0 = 1.4
  l_max_its = 100
  eig_check_tol = 1e-7

  solve_type = 'NEWTON'
  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_linesearch_monitor'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       NONZERO               superlu_dist'

## Use the settings below instead if running on a desktop/small cluster
#  petsc_options_iname = '-pc_type -sub_pc_type -ksp_gmres_restart -pc_asm_overlap -sub_pc_factor_shift_type'
#  petsc_options_value = 'asm      lu           200                1               NONZERO'

  line_search = none
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Postprocessors]
  [./bnorm]
    type = ElmIntegTotFissNtsPostprocessor
    execute_on = linear
  [../]
  [./tot_fissions]
    type = ElmIntegTotFissPostprocessor
    execute_on = linear
  [../]
  [./powernorm]
    type = ElmIntegTotFissHeatPostprocessor
    execute_on = linear
  [../]
  [./group1norm]
    type = ElementIntegralVariablePostprocessor
    variable = group1
    execute_on = linear
  [../]
  [./group1max]
    type = NodalMaxValue
    variable = group1
    execute_on = timestep_end
  [../]
  [./group1diff]
    type = ElementL2Diff
    variable = group1
    execute_on = 'linear timestep_end'
    use_displaced_mesh = false
  [../]
  [./group2norm]
    type = ElementIntegralVariablePostprocessor
    variable = group2
    execute_on = linear
  [../]
  [./group2max]
    type = NodalMaxValue
    variable = group2
    execute_on = timestep_end
  [../]
  [./group2diff]
    type = ElementL2Diff
    variable = group2
    execute_on = 'linear timestep_end'
    use_displaced_mesh = false
  [../]
[]

[Outputs]
  perf_graph = true
  print_linear_residuals = true
  [./exodus]
    type = Exodus
  [../]
[]

[Debug]
  show_var_residual_norms = true
[]

The radial grid looks like this:

But it fails to converge, the log is as follows:

(moose) jixu@DESKTOP-9HFBUBM:~/projects/moltres/lfr/cnrs_midified/phase1$ mpiexec -n 2 ~/projects/moltres/moltres-opt -i eigen.i


*** Warning, This code is deprecated and will be removed in future versions:
***** Invalid Object: NodalMaxValue *****
Expired on Mon Sep 27 01:00:00 2021
Update your application using the 'NodalExtremeValue' object



*** Info ***
'InversePowerMethod' executioner is deprecated in favor of 'Eigenvalue' executioner.
Few parameters such as 'bx_norm', 'k0', 'xdiff', 'max_power_iterations', 'min_power_iterations', 'eig_check_tol', 'sol_check_tol', and 'output_before_normalization' are no longer supported.
However, 'Eigenvalue' executioner supports more solving options by interfacing SLEPc.



*** Warning, This code is deprecated and will be removed in future versions:
AddVariableAction::determineType() is deprecated. Use AddVariableAction::variableType() instead.


  Initializing
    Initializing Equation Systems                                                        [  5.98 s] [  830 MB]
  Finished Initializing                                                                  [  5.99 s] [  831 MB]
Finished Setting Up                                                                      [  6.40 s] [  860 MB]
Framework Information:
MOOSE Version:           git commit 74571c498f on 2022-06-16
LibMesh Version:
PETSc Version:           3.16.5
SLEPc Version:           3.16.2
Current Time:            Tue Aug  2 16:28:21 2022
Executable Timestamp:    Fri Jul  8 14:37:35 2022

Parallelism:
  Num Processors:          2
  Num Threads:             1

Mesh:
  Parallel Type:           replicated
  Mesh Dimension:          3
  Spatial Dimension:       3
  Nodes:
    Total:                 16149
    Local:                 8421
    Min/Max/Avg:           7728/8421/8074
  Elems:
    Total:                 14880
    Local:                 7440
    Min/Max/Avg:           7440/7440/7440
  Num Subdomains:          3
  Num Partitions:          2
  Partitioner:             metis

Nonlinear System:
  Num DOFs:                186174
  Num Local DOFs:          95166
  Variables:               { "group1" "group2" "group3" "group4" "group5" "group6" } { "pre1" "pre2"
                             "pre3" "pre4" "pre5" "pre6" }
  Finite Element Types:    "LAGRANGE" "MONOMIAL"
  Approximation Orders:    "FIRST" "CONSTANT"

Execution Information:
  Executioner:             InversePowerMethod
  Solver Mode:             NEWTON
  MOOSE Preconditioner:    SMP

LEGACY MODES ENABLED:
 This application uses the legacy material output option: material properties are output only on TIMESTEP_END, not INITIAL. To remove this message, set 'use_legacy_material_output' to false in this application. If there are gold output files that contain material property output for which output occurs on INITIAL, then these will generate diffs due to zero values being stored, and these tests should be re-golded.


|Bx| = 1.4000000000

Time Step 0, time = 0

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| time           | bnorm          | group1diff     | group1max      | group1norm     | group2diff     | group2max      | group2norm     | powernorm      | tot_fissions   |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   1.400000e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   1.538358e-11 |   4.833036e-01 |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+


Reporter Values:
+----------------+------------------------+
| time           | Executioner/eigenvalue |
+----------------+------------------------+
|   0.000000e+00 |           1.400000e+00 |
+----------------+------------------------+



*** Warning ***
/home/jixu/projects/moltres/lfr/cnrs_midified/phase1/eigen.i:24: unused parameter 'Nt/transient'


 Power iterations starts
 ________________________________________________________________________________
 Power iteration= 0
 0 Nonlinear |R| = 1.229227e+01
    |residual|_2 of individual variables:
                     group1: 4.43834
                     group2: 6.68472
                     group3: 2.86998
                     group4: 3.25508
                     group5: 3.70077
                     group6: 7.36115
                     pre1:   0.00161329
                     pre2:   0.0121789
                     pre3:   0.00996549
                     pre4:   0.0230553
                     pre5:   0.0134822
                     pre6:   0.00484308

    Computing Jacobian                                                                   [  6.36 s] [  340 MB]
      0 Linear |R| = 1.229227e+01
      1 Linear |R| = 1.229227e+01
      2 Linear |R| = 1.229227e+01
      3 Linear |R| = 1.229227e+01
      4 Linear |R| = 1.229226e+01
      5 Linear |R| = 1.229226e+01
      6 Linear |R| = 1.229226e+01
      7 Linear |R| = 1.229226e+01
      8 Linear |R| = 1.229226e+01
      9 Linear |R| = 1.229225e+01
     10 Linear |R| = 1.229225e+01
     11 Linear |R| = 1.229225e+01
     12 Linear |R| = 1.229225e+01
     13 Linear |R| = 1.229225e+01
     14 Linear |R| = 1.229224e+01
     15 Linear |R| = 1.229224e+01
     16 Linear |R| = 1.229224e+01
     17 Linear |R| = 1.229223e+01
     18 Linear |R| = 1.229223e+01
     19 Linear |R| = 1.229223e+01
     20 Linear |R| = 1.229223e+01
     21 Linear |R| = 1.229223e+01
     22 Linear |R| = 1.229222e+01
     23 Linear |R| = 1.229222e+01
     24 Linear |R| = 1.229222e+01
     25 Linear |R| = 1.229222e+01
     26 Linear |R| = 1.229222e+01
     27 Linear |R| = 1.229222e+01
     28 Linear |R| = 1.229221e+01
     29 Linear |R| = 1.229220e+01
     30 Linear |R| = 1.229220e+01
     31 Linear |R| = 1.229220e+01
     32 Linear |R| = 1.229220e+01
     33 Linear |R| = 1.229220e+01
     34 Linear |R| = 1.229220e+01
     35 Linear |R| = 1.229220e+01
     36 Linear |R| = 1.229220e+01
     37 Linear |R| = 1.229220e+01
     38 Linear |R| = 1.229220e+01
     39 Linear |R| = 1.229220e+01
     40 Linear |R| = 1.229220e+01
     41 Linear |R| = 1.229220e+01
     42 Linear |R| = 1.229220e+01
     43 Linear |R| = 1.229220e+01
     44 Linear |R| = 1.229220e+01
     45 Linear |R| = 1.229220e+01
     46 Linear |R| = 1.229220e+01
     47 Linear |R| = 1.229220e+01
     48 Linear |R| = 1.229220e+01
     49 Linear |R| = 1.229220e+01
     50 Linear |R| = 1.229220e+01
     51 Linear |R| = 1.229220e+01
     52 Linear |R| = 1.229220e+01
     53 Linear |R| = 1.229220e+01
     54 Linear |R| = 1.229220e+01
     55 Linear |R| = 1.229220e+01
     56 Linear |R| = 1.229220e+01
     57 Linear |R| = 1.229220e+01
     58 Linear |R| = 1.229220e+01
     59 Linear |R| = 1.229220e+01
     60 Linear |R| = 1.229220e+01
     61 Linear |R| = 1.229220e+01
     62 Linear |R| = 1.229220e+01
     63 Linear |R| = 1.229220e+01
     64 Linear |R| = 1.229220e+01
     65 Linear |R| = 1.229220e+01
     66 Linear |R| = 1.229220e+01
     67 Linear |R| = 1.229220e+01
     68 Linear |R| = 1.229220e+01
     69 Linear |R| = 1.229220e+01
     70 Linear |R| = 1.229220e+01
     71 Linear |R| = 1.229220e+01
     72 Linear |R| = 1.229220e+01
     73 Linear |R| = 1.229220e+01
     74 Linear |R| = 1.229220e+01
     75 Linear |R| = 1.229220e+01
     76 Linear |R| = 1.229220e+01
     77 Linear |R| = 1.229220e+01
     78 Linear |R| = 1.229220e+01
     79 Linear |R| = 1.229220e+01
     80 Linear |R| = 1.229220e+01
     81 Linear |R| = 1.229220e+01
     82 Linear |R| = 1.229220e+01
     83 Linear |R| = 1.229220e+01
     84 Linear |R| = 1.229220e+01
     85 Linear |R| = 1.229220e+01
     86 Linear |R| = 1.229220e+01
     87 Linear |R| = 1.229220e+01
     88 Linear |R| = 1.229220e+01
     89 Linear |R| = 1.229220e+01
     90 Linear |R| = 1.229220e+01
     91 Linear |R| = 1.229220e+01
     92 Linear |R| = 1.229220e+01
     93 Linear |R| = 1.229220e+01
     94 Linear |R| = 1.229220e+01
     95 Linear |R| = 1.229220e+01
     96 Linear |R| = 1.229220e+01
     97 Linear |R| = 1.229220e+01
     98 Linear |R| = 1.229220e+01
     99 Linear |R| = 1.229220e+01
    100 Linear |R| = 1.229220e+01
  Linear solve did not converge due to DIVERGED_ITS iterations 100
 1 Nonlinear |R| = 1.229220e+01
    |residual|_2 of individual variables:
                     group1: 4.43824
                     group2: 6.68466
                     group3: 2.86998
                     group4: 3.25508
                     group5: 3.70078
                     group6: 7.36115
                     pre1:   0.00161329
                     pre2:   0.0121789
                     pre3:   0.00996549
                     pre4:   0.0230553
                     pre5:   0.0134822
                     pre6:   0.00484308
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
  Finished Solving                                                                       [189.43 s] [ 5314 MB]

 +================+=====================+=====================+
 | iteration      | eigenvalue          | solution_difference |
 +================+=====================+=====================+
 |              0 |      1.40000000e+00 |      8.75503036e-08 |
 +================+=====================+=====================+
 Chebyshev step: 0
 ________________________________________________________________________________

*******************************************************
 Eigenvalue = 1.3999999986
*******************************************************

Outlier Variable Residual Norms:
  group2: 6.684662e+00
  group6: 7.361151e+00

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| time           | bnorm          | group1diff     | group1max      | group1norm     | group2diff     | group2max      | group2norm     | powernorm      | tot_fissions   |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   1.400000e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   1.538358e-11 |   4.833036e-01 |
|   1.000000e+00 |   1.400000e+00 |   8.755030e-08 |   7.209564e-02 |   8.432560e+00 |   9.110021e-08 |   7.209589e-02 |   8.432560e+00 |   1.538358e-11 |   4.833036e-01 |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+


Reporter Values:
+----------------+------------------------+
| time           | Executioner/eigenvalue |
+----------------+------------------------+
|   0.000000e+00 |           1.400000e+00 |
|   1.000000e+00 |           1.400000e+00 |
+----------------+------------------------+

 Solution is rescaled with factor 1.38004e+15 for normalization!

Outlier Variable Residual Norms:
  group2: 6.684662e+00
  group6: 7.361151e+00

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
| time           | bnorm          | group1diff     | group1max      | group1norm     | group2diff     | group2max      | group2norm     | powernorm      | tot_fissions   |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   1.400000e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   7.797104e-01 |   0.000000e+00 |   8.432560e+00 |   1.538358e-11 |   4.833036e-01 |
|   1.000000e+00 |   1.400000e+00 |   8.755030e-08 |   7.209564e-02 |   8.432560e+00 |   9.110021e-08 |   7.209589e-02 |   8.432560e+00 |   1.538358e-11 |   4.833036e-01 |
|   2.000000e+00 |   1.932061e+15 |   1.076034e+15 |   9.949510e+13 |   1.163730e+16 |   1.076034e+15 |   9.949544e+13 |   1.163730e+16 |   2.123000e+04 |   6.669799e+14 |
+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+----------------+


Reporter Values:
+----------------+------------------------+
| time           | Executioner/eigenvalue |
+----------------+------------------------+
|   0.000000e+00 |           1.400000e+00 |
|   1.000000e+00 |           1.400000e+00 |
|   2.000000e+00 |           1.400000e+00 |
+----------------+------------------------+



Performance Graph:
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------
|                                Section                               | Calls |   Self(s)  |   Avg(s)   |    %   | Mem(MB) |  Total(s)  |   Avg(s)   |    %   | Mem(MB) |
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------
| MoltresApp (main)                                                    |     1 |      0.010 |      0.010 |   0.00 |       3 |    198.585 |    198.585 | 100.00 |    6193 |
|   Action::SetupMeshAction::Mesh::SetupMeshAction::act::setup_mesh    |     1 |      0.001 |      0.001 |   0.00 |       0 |      0.001 |      0.001 |   0.00 |       0 |
|   Action::SetupMeshAction::Mesh::SetupMeshAction::act::set_mesh_base |     2 |      0.001 |      0.001 |   0.00 |       0 |      0.167 |      0.083 |   0.08 |      14 |
|   FEProblem::computeUserObjects                                      |     5 |      1.056 |      0.211 |   0.53 |       2 |      1.056 |      0.211 |   0.53 |       2 |
|   FEProblem::outputStep                                              |     3 |      0.001 |      0.000 |   0.00 |       0 |      0.506 |      0.169 |   0.25 |       6 |
|     Exodus::outputStep                                               |     3 |      0.460 |      0.153 |   0.23 |       6 |      0.460 |      0.153 |   0.23 |       6 |
|   FEProblem::solve                                                   |     1 |    175.603 |    175.603 |  88.43 |    4956 |    189.432 |    189.432 |  95.39 |    5314 |
|     FEProblem::computeResidualInternal                               |     1 |      0.000 |      0.000 |   0.00 |       0 |      2.460 |      2.460 |   1.24 |       1 |
|       FEProblem::computeUserObjects                                  |     1 |      0.335 |      0.335 |   0.17 |       0 |      0.335 |      0.335 |   0.17 |       0 |
|     FEProblem::computeResidualInternal                               |     2 |      0.002 |      0.001 |   0.00 |       0 |      4.943 |      2.471 |   2.49 |      14 |
|       FEProblem::computeUserObjects                                  |     2 |      0.692 |      0.346 |   0.35 |       6 |      0.692 |      0.346 |   0.35 |       6 |
|     FEProblem::computeJacobianInternal                               |     1 |      0.000 |      0.000 |   0.00 |       0 |      6.361 |      6.361 |   3.20 |     340 |
|   InversePowerMethod::final                                          |     1 |      0.000 |      0.000 |   0.00 |       0 |      0.001 |      0.001 |   0.00 |       0 |
|     FEProblem::outputStep                                            |     1 |      0.000 |      0.000 |   0.00 |       0 |      0.001 |      0.001 |   0.00 |       0 |
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Finished Executing                                                        

               [192.17 s] [ 5330 MB]

I want to know how to deal with it. Any suggestions are sincerely appreciated and all files are attached.
files.zip

[smpark7]

My computer doesn't have enough memory to run the 3D simulation so I can't comment on your 3D mesh, but the mesh looks alright. I found two minor errors in your input file according to your description of what you intend to achieve:

1. You have to define pre_blocks = fuel in the Nt section because you have precursors in the fuel block only.
2. Similarly, you have to define block = fuel in the Precursors section.

Here's a working example with a 2D mesh that I created with MeshGenerator from your 3D mesh dimensions: eigen.txt
The 2D example has no vacuum boundaries because it doesn't contain the top and bottom surfaces. I got k=1.561946, which would probably drop to 1.3-1.4 when you run the 3D model instead. I hope this helps with your simulations!

[js-jixu]

Thank you for your helpful advice, Sun Myung. Your are so experienced that you can guess that the eigenvalue of 3D model is between 1.3 and 1.4. Indeed, the eigenvalue calculated by OpenMC is 1.34421 +/- 0.00022, while Moltres's results is 1.334104. Their difference is about pcm, which is a small error for a 3D model.

After I changed the -pc_type parameter to asm yesterday, it worked successfully. The results obtained are very close to those of OpenMC.

So I thought maybe not defining pre_blocks=fuel in Nt section and defining block=fuel in Precursors section would also get the correct result. Because although I have defined precursors in cool, clad and fuel (I know that precursors are only in fuel), the group constants BETA_EFF, CHI_D, CHI_P, CHI_T, FISS, KAPPA, LAMBDA, NSF of clad and cool are all equal to 0, So the precursors of the clad and cool areas should be equal to 0. Do you think there is a problem with this understanding?

After I defined pre_blocks = fuel in Nt section and block = fuel in Precursors section, the results of Moltres also is 1.334104. When they are not defined, it takes 150s to calculate one step of nonlinear iteration, and it takes 100s to calculate one step after defining.

When I plot a line from (0, -0.57, 45) to (0, 0.57, 45), the results shows that the neutron flux and precursors are uniformly distributed in the radial direction and they only vary in the axial direction. The distribution situation is not that the neutron flux in the fuel region is relatively high, and the neutron flux in the other places is relatively small. Do you think this result is reasonable?

There are new input file, mesh files and output file.
new files.zip

[smpark7]

The flux variation in the radial direction is very small compared to the axial direction because your pin is only 1.14cm across and 45 cm long. Here are the group 1 and 6 fluxes from the 2D model I ran with your cross sections.

You should compare them with your OpenMC results to check if they're correct.

Also, you can initialize precursors in the clad and coolant region, but that will increase your problem size because you're solving for more values. This is why it takes longer per iteration if you have precursors in those regions. Moltres will try to solve for their values because Moltres doesn't know that they're always zero.

[js-jixu]

Yeah, you're right, Park. If I want to do a coupling case similar to \moltres\problems\2021-cnrs-benchmark\phase-1\full-coupling.i, how could I export the EigenValue data of the neutron field calculated at each step. In another word, I want to get the EigenValue for each time step in a table, how can I do it?