Why I use automatic_scaling, Convergence is worse?

[wangzhaohao]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello, MOOSE developers.
I am going to use SUPG Method in a diffusion-advection problem. when I code this, follwing is the PDE.
$$\rho C \frac{\partial T}{\partial t} - \nabla \cdot (k \nabla T) - q = 0$$
$$\frac{\partial p}{\partial t} + \nabla \cdot [(1 - p) p \vec{v} - \nu \nabla p] = 0$$
$$0 = \rho \cdot \Delta$$
So my Input file is

[GlobalParams]
  coord_type = RZ
  displacements = 'disp_x'
[]

[Mesh]
  [square]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 50
    xmax = 2.675e-3
  []
[]

[Variables]
  [disp_x]
  []
  [scalar_strain_yy]
    order = FIRST
    family = SCALAR
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300
  []
  [porosity]
    initial_condition = 0.15
  []
[]

[AuxVariables]
  [disp_z]
  []
  [velocity]
    family = MONOMIAL_VEC
    order = FIRST
  []
  [grad_temp_x]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [velocity]
    type = PoreVelocityAux
    variable = velocity
    temperature = temp
    plutonium_concentration = 0
  []
  [grad_temp_x]
    type = VariableGradientComponent
    variable = grad_temp_x
    component = x
    gradient_variable = temp
  []
[]

[Postprocessors]
  [react_z]
    type = MaterialTensorIntegral
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
  []
[]

[Kernels]
  inactive = 'heatsourceconst'
  [heatconduction]
    type = HeatConduction
    variable = temp
  []
  [heatsourceconst]
    type = HeatSource
    function = source_fun
    variable = temp
  []
  [heatsource]
    type = OxidalHeatSource
    function = source_fun
    porosity = porosity
    porosity_initial = 0.15
    variable = temp
  []
  [heatconductiontimederivative]
    type = HeatConductionTimeDerivative
    specific_heat = 120
    variable = temp
  []

  [porositymigration]
    type = PoreMigration
    variable = porosity
    temperature = temp
    plutonium_concentration = 0
  []
  [porosityadvectionSUPG]
    type = PoreAdvectionSUPG
    temperature = temp
    variable = porosity
    plutonium_concentration = 0
  []
  [poretimederivative]
    type = TimeDerivative
    variable = porosity
  []
  [porositytimederivativeSUPG]
    type = PoreTimeDerivativeSUPG
    temperature = temp
    variable = porosity
  []
[]

[Physics]

  [SolidMechanics]

    [QuasiStatic]
      [all]
        strain = SMALL
        incremental = true
        generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy strain_zz'
        planar_formulation = GENERALIZED_PLANE_STRAIN
        eigenstrain_names = eigenstrain
        scalar_out_of_plane_strain = scalar_strain_yy
        temperature = temp
        out_of_plane_direction = z
      []
    []
  []
[]

[Functions]
  [temp_bc]
    type = PiecewiseLinear
    xy_data = '0 0.63e3
              1e4 1.3e3'
  []
  [pressure_bc]
    type = PiecewiseLinear
    xy_data = '0 0.2e6
              1e4 0.4e6'
  []
  [source_fun]
    type = PiecewiseLinear
    xy_data = '0 0
              1e4 2.2253216e9' #5e4W/m
  []
[]

[BCs]
  [bottomx]
    type = DirichletBC
    boundary = 0
    variable = disp_x
    value = 0.0
  []
  [temp]
    type = FunctionDirichletBC
    variable = temp
    boundary = 1
    function = temp_bc
  []
  [outpressure]
    type = Pressure
    boundary = 1
    variable = disp_x
    function = pressure_bc
  []
[]

[Materials]
  [elastic_tensor]
    type = ComputeIsotropicElasticityTensor
    poissons_ratio = 0.345
    youngs_modulus = 200e9
  []
  [thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    temperature = temp
    thermal_expansion_coeff = 1e-6
    stress_free_temperature = 300
    eigenstrain_name = eigenstrain
  []
  [stress]
    type = ComputeStrainIncrementBasedStress
  []
  [density]
    type = Density
    density = 1.0662e4
  []
  [mamox]
    type = MAMOXThermal
    temperature = temp
    porosity = porosity
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  line_search = none
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  automatic_scaling = true
  scaling_group_variables = 'disp_x; temp; porosity'

  # controls for linear iterations
  l_max_its = 100
  l_tol = 1e-4

  # controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-5

  # time control
  start_time = 0.0
  dt = 40
  dtmin = 1.0
  end_time = 1e4
  num_steps = 5000
[]

[Outputs]
  [exodus]
    type = Exodus
    output_material_properties = true
    show_material_properties = 'thermal_conductivity'
  []
[]

If I don't use the automatic_scaling = true scaling_group_variables = 'disp_x; temp; porosity', the residual in the last step is like this

Time Step 251, time = 10000, dt = 20
 0 Nonlinear |R| = 1.028307e+04
      0 Linear |R| = 1.028307e+04
      1 Linear |R| = 8.921250e+02
      2 Linear |R| = 2.291096e+01
      3 Linear |R| = 1.745521e+01
      4 Linear |R| = 1.601492e+01
      5 Linear |R| = 9.808384e+00
      6 Linear |R| = 9.610114e+00
      7 Linear |R| = 6.457256e+00
      8 Linear |R| = 5.269686e+00
      9 Linear |R| = 3.571940e+00
     10 Linear |R| = 3.205730e+00
     11 Linear |R| = 2.177902e+00
     12 Linear |R| = 2.091932e+00
     13 Linear |R| = 1.138776e+00
     14 Linear |R| = 1.097427e+00
     15 Linear |R| = 8.563212e-01
 1 Nonlinear |R| = 9.580684e+00
      0 Linear |R| = 9.580684e+00
      1 Linear |R| = 9.301138e+00
      2 Linear |R| = 9.300915e+00
      3 Linear |R| = 9.170254e+00
      4 Linear |R| = 7.566562e+00
      5 Linear |R| = 7.185175e+00
      6 Linear |R| = 4.388656e+00
      7 Linear |R| = 4.386418e+00
      8 Linear |R| = 2.456251e+00
      9 Linear |R| = 1.690797e+00
     10 Linear |R| = 1.272286e+00
     11 Linear |R| = 1.221710e+00
     12 Linear |R| = 1.213552e+00
     13 Linear |R| = 4.786498e-01
     14 Linear |R| = 4.148035e-01
     15 Linear |R| = 2.105089e-01
     16 Linear |R| = 2.104989e-01
     17 Linear |R| = 1.100865e-01
     18 Linear |R| = 8.341243e-02
     19 Linear |R| = 4.653824e-02
     20 Linear |R| = 2.011911e-02
     21 Linear |R| = 1.194144e-02
     22 Linear |R| = 7.006019e-03
     23 Linear |R| = 6.637877e-03
     24 Linear |R| = 4.550686e-03
     25 Linear |R| = 2.480034e-03
     26 Linear |R| = 1.546153e-03
     27 Linear |R| = 7.928538e-04
 2 Nonlinear |R| = 1.277763e+00
      0 Linear |R| = 1.277763e+00
      1 Linear |R| = 1.187731e+00
      2 Linear |R| = 1.169138e+00
      3 Linear |R| = 8.000056e-01
      4 Linear |R| = 5.099548e-01
      5 Linear |R| = 5.094007e-01
      6 Linear |R| = 2.483796e-01
      7 Linear |R| = 2.460117e-01
      8 Linear |R| = 1.194662e-01
      9 Linear |R| = 1.133298e-01
     10 Linear |R| = 4.650997e-02
     11 Linear |R| = 3.583607e-02
     12 Linear |R| = 1.449329e-02
     13 Linear |R| = 1.068608e-02
     14 Linear |R| = 6.163695e-03
     15 Linear |R| = 6.135168e-03
     16 Linear |R| = 2.384477e-03
     17 Linear |R| = 2.370099e-03
     18 Linear |R| = 1.001684e-03
     19 Linear |R| = 9.856626e-04
     20 Linear |R| = 4.689908e-04
     21 Linear |R| = 3.381278e-04
     22 Linear |R| = 1.462069e-04
     23 Linear |R| = 6.179482e-05
 3 Nonlinear |R| = 6.335172e-02
 Solve Converged!

Outlier Variable Residual Norms:
  disp_x: 6.335137e-02

but when i used

Framework Information:
MOOSE Version:           git commit 2e75e9d1e3 on 2024-11-26
LibMesh Version:         
PETSc Version:           3.22.1
SLEPc Version:           3.22.1
Current Time:            Mon Dec  2 11:27:08 2024
Executable Timestamp:    Sun Dec  1 17:13:12 2024

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:          1
  Spatial Dimension:       1
  Nodes:                   51
  Elems:                   50
  Num Subdomains:          1

Nonlinear System:
  Num DOFs:                154
  Num Local DOFs:          154
  Variables:               "disp_x" "scalar_strain_yy" { "temp" "porosity" } 
  Finite Element Types:    "LAGRANGE" "SCALAR" "LAGRANGE" 
  Approximation Orders:    "FIRST" "FIRST" "FIRST" 

Auxiliary System:
  Num DOFs:                651
  Num Local DOFs:          651
  Variables:               "disp_z" "velocity" { "grad_temp_x" "stress_xx" "stress_xy" "stress_yy" "stress_zz" 
                             "strain_xx" "strain_xy" "strain_yy" "strain_zz" "thermal_conductivity" } 
  Finite Element Types:    "LAGRANGE" "MONOMIAL_VEC" "MONOMIAL" 
  Approximation Orders:    "FIRST" "FIRST" "CONSTANT" 

Execution Information:
  Executioner:             Transient
  TimeStepper:             ConstantDT
  TimeIntegrator(s):       ImplicitEuler
  Solver Mode:             Preconditioned JFNK
  PETSc Preconditioner:    lu 
  MOOSE Preconditioner:    SMP


Time Step 0, time = 0

Postprocessor Values:
+----------------+----------------+
| time           | react_z        |
+----------------+----------------+
|   0.000000e+00 |   0.000000e+00 |
+----------------+----------------+


Scalar Variable Values:
+----------------+------------------+
| time           | scalar_strain_yy |
+----------------+------------------+
|   0.000000e+00 |     0.000000e+00 |
+----------------+------------------+


Time Step 1, time = 40, dt = 40

Performing automatic scaling calculation

 0 Nonlinear |R| = 6.263956e+09
    Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
                   PC failed due to FACTOR_NUMERIC_ZEROPIVOT 


*** Warning ***
/home/wang/projects/miracle/test/tests/kernels/restructure/restructure_SUPG.i:194.3:
The following warning occurred in the MaterialBase 'mamox' of type MAMOXThermal.

WARNING: In MAMOXThermal: negative temperature!
         Reseting to zero.
        _qp: 0
        temp: -inf
        elem: 0
        proc: 0

I also try to use analyzejacobian to check the Jacobian.

** PETSc DEPRECATION WARNING ** : the option -snes_test_jacobian_display_threshold is deprecated as of version 3.13 and will be removed in a future release.
   Use the option -snes_test_jacobian instead. -snes_test_jacobian accepts an optional threshold (since v3.10) (Silence this warning with -options_suppress_deprecated_warnings)
  ---------- Testing Jacobian -------------
  Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is
    O(1.e-8), the hand-coded Jacobian is probably correct.
  ||J - Jfd||_F/||J||_F = 3.52237e-08, ||J - Jfd||_F = 2.77796e+08
  Hand-coded Jacobian ----------
Mat Object: (nl0_) 1 MPI process

So how to fix this problem, thanks for you help in advance.
Zhaohao Wang.

[GiudGiud]

hello

you can try

  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'

to get around a zero-pivot in LU

This is fairly loose

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-5

so you might want to tighten it or be mindful if it causes trouble later on

Guillaume

[wangzhaohao]

Thanks for your quick reply, I have tried the method you provided, and reduse the nl_rel_tol nl_abs_tol
the result seems not well.
When I use the automatic_scaling, It's over very quickly, like noting happened. The variable 'porosity' doesn't have any redistribution. and the variable 'temperature' is quite strange.

If close the automatic_scaling,

48 Nonlinear |R| = 8.139512e-04
      0 Linear |R| = 8.139512e-04
      1 Linear |R| = 8.133132e-04
      2 Linear |R| = 8.118787e-04
      3 Linear |R| = 7.965513e-04
      4 Linear |R| = 7.382100e-04
      5 Linear |R| = 7.006991e-04
      6 Linear |R| = 2.310533e-04
      7 Linear |R| = 1.289893e-04
      8 Linear |R| = 3.198177e-05
      9 Linear |R| = 1.427464e-05
     10 Linear |R| = 4.980067e-06
     11 Linear |R| = 1.327768e-06
     12 Linear |R| = 6.281298e-07
     13 Linear |R| = 2.168976e-07
     14 Linear |R| = 7.426662e-08
    |residual|_2 of individual variables:
           disp_x:           0.000828241
           scalar_strain_yy: 7.27083e-09
           temp:             3.24275e-06
           porosity:         6.98804e-12
49 Nonlinear |R| = 8.282470e-04
      0 Linear |R| = 8.282470e-04
      1 Linear |R| = 8.264407e-04
      2 Linear |R| = 8.238493e-04
      3 Linear |R| = 8.086120e-04
      4 Linear |R| = 7.590600e-04
      5 Linear |R| = 6.051407e-04
      6 Linear |R| = 2.515633e-04
      7 Linear |R| = 1.100417e-04
      8 Linear |R| = 2.456599e-05
      9 Linear |R| = 1.323413e-05
     10 Linear |R| = 3.790405e-06
     11 Linear |R| = 1.344224e-06
     12 Linear |R| = 6.340606e-07
     13 Linear |R| = 2.036445e-07
     14 Linear |R| = 7.245698e-08
    |residual|_2 of individual variables:
           disp_x:           0.000923801
           scalar_strain_yy: 5.78936e-10
           temp:             3.43356e-06
           porosity:         4.3539e-12
50 Nonlinear |R| = 9.238078e-04
  Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 50
 Solve Did NOT Converge!
Aborting as solve did not converge


*** ERROR ***
The following error occurred in the TimeStepper 'ConstantDT' of type ConstantDT.

Solve failed and timestep already at or below dtmin, cannot continue!

[GiudGiud]

If you use AD (by using an ADKernel and consistently using AD numbers when computing the residual) you don’t need to write any of the computerJacobian routines.

[wangzhaohao]

Hello, I find a way how to fix this problem in Converge. May it's a diffucult in using FEM to solve fluid problems.
Look back the PDE. The origin PDE about pore redistribution model is $\frac{dp}{dt}+\Delta \cdot (\bar{v}p -\mu \Delta p) = 0$, it's hard to solve with FEM, particulary $\bar{v}$ change with temperature and gradient. So we modify this PDE to $\frac{dp}{dt}+\Delta \cdot (\bar{v}p(1-p) -\mu \Delta p) = 0$, compared with origin, which has batter converage.
About the automatic_scaling, I find that in this problem it led MOOSE in the wrong direction. So I close it, and the converage is convergence is barely satisfactory. Hopefully it will converge in multiple dimensions.

    |residual|_2 of individual variables:
           disp_x:           0.976429
           scalar_strain_yy: 0.000333498
           temp:             0.00533607
           porosity:         2.01842e-14
 2 Nonlinear |R| = 9.764435e-01
      0 Linear |R| = 9.764435e-01
      1 Linear |R| = 2.234593e-02
      2 Linear |R| = 5.979464e-03
      3 Linear |R| = 3.565020e-04
      4 Linear |R| = 2.572591e-04
      5 Linear |R| = 5.715189e-05
    |residual|_2 of individual variables:
           disp_x:           7.62918e-05
           scalar_strain_yy: 2.44278e-05
           temp:             5.11907e-05
           porosity:         2.74437e-16
 3 Nonlinear |R| = 9.506656e-05
 Solve Converged!

I have not studied fluid mechanics systematically, and the above are just my personal opinions. If you have other ideas, please let me know, thank you very much.

[GiudGiud]

you could try using multi-system with picard iterations. So set the variables to be solved in different numerical systems (solver_sys parameter in each Variable),add the list of systems in nl_sys_names in Problem, and set additional parameters to the executioner to trigger the fixed point iteration

this will replace Newton on all 4 variables with Newton on a reduced set of variables and fixed point iterations.
If there are variables which are tightly coupled, I recommend they be kept in the same nonlinear system

[wangzhaohao]

Thanks, it looks different from the MultiApp I've used. I've already gone through this part, and I'll give it a try later.

[GiudGiud]

It's fairly similar in terms of convergence properties, except you can use relaxation with the multiapps and it is not implemented yet with the multi-system.
The performance of multi-system should be a tad better, considering you do not have to transfer variables.