Maybe there are some bugs in navier stokes modules

[suqingji]

Dear developers,

I am developing the RANS k-omega turbulence model based on MOOSE FEM. Now my k-omega codes can be used to simulate some simple cases. But in other complex cases, the convergence of the code is difficult. So I find some bugs(maybe which is a normal phenomenon but I think it is abnormal) in laminar simulation when I debug my turbulence codes.

Here is the list about the bugs(Please allow me to call it a bug).

1. In a 3D-lid-driven test case, I found that the convergence progress is different on different scale grid unless the variables are scaled.

for example,
If I change the L params by multiplying factor 1000 with type = TransformGenerator, at the same time I change the \mu(Dynamic viscosity) params by multiplying factor 1000, so the Reynolds number and Euler number are the same. And I also make time step dt becoming 1000 times. So the Strauhal number is the same, too. In other words, the above two problems are the same. So that the convergence progress should be the same in theory. But It is amazing to me that the convergence progress is different.

Maybe you will tell me that it does not matter. But in fact in complex problems especially in Thermal Hydraulics, the problem may become Solve Did NOT Converge!

tips: In the above problem, I use AD rather than by hand to compute Jacobian.

2. Also in a 3D-lid-driven test case

Now I use by hand rather than AD to compute Jacobian. And the conditions are exactly the same as list 1. So that the solution should be the same in theory. And the converge or not converge should be the same in theory. But It is amazing to me that the convergence is different. In AD the problem is converged but in by hand the problem is not converged. If I change the SUPG/PSPG params alpha in the input file the problem becomes converged. Those are unbelievable.

It sounds crazy, and it doesn't make sense in theory, but I've actually found these problems, so it feels like some bugs to me.

Best regards.

QINGJI SU

[lindsayad]

I would have to see the hand-coded and AD inputs in order to understand why you are getting different behavior between the two cases.

Regarding the scaling ... I believe what you are seeing there is a consequence of the fact that degrees of freedom controlled by Dirichlet boundary conditions are included in the nonlinear system. So as you scale your physics rows, your Dirichlet rows remain unscaled and so you have different operators

[suqingji]

About the second problem, Here are my test case input files.
AD:

const_num = 1.e0
const_dt = 1.e-1
const_mu = ${const_num}
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    xmax = ${const_num}
    ymin = 0
    ymax = ${const_num}
    zmin = 0
    zmax = ${const_num}
    nx = 10
    ny = 10
    nz = 10
  []
  [./corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = gen
  [../]
[]
[Variables]
  [./velocity]
    family = LAGRANGE_VEC
  [../]
  [./p]
  [../]
[]
[ICs]
  [velocity]
    type = VectorConstantIC
    x_value = 1e-15
    y_value = 1e-15
    variable = velocity
  []
[]
[Kernels]
  [./mass]
    type = INSADMass
    variable = p
  [../]
  [./mass_pspg]
    type = INSADMassPSPG
    variable = p
  [../]
  [./momentum_time]
    type = INSADMomentumTimeDerivative
    variable = velocity
  [../]
  [./momentum_convection]
    type = INSADMomentumAdvection
    variable = velocity
  [../]
  [./momentum_viscous]
    type = INSADMomentumViscous
    variable = velocity
  [../]
  [./momentum_pressure]
    type = INSADMomentumPressure
    variable = velocity
    pressure = p
    integrate_p_by_parts = true
  [../]
  [./momentum_supg]
    type = INSADMomentumSUPG
    variable = velocity
    velocity = velocity
  [../]
[]
[BCs]
  [./no_slip]
    type = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'bottom right left front back'
  [../]
  [./lid]
    type = VectorFunctionDirichletBC
    variable = velocity
    boundary = 'top'
    function_x = '1'
  [../]
  [./pressure_pin]
    type = DirichletBC
    variable = p
    boundary = 'pinned_node'
    value = 0
  [../]
[]
[Materials]
  [./const]
    type = ADGenericConstantMaterial
    prop_names = 'rho      mu                 cp k'
    prop_values = '1e4     ${const_mu}       1  1e4'
  [../]
  [ins_mat]
    type = INSADTauMaterial
    velocity = velocity
    pressure = p
    alpha = 1
  []
[]
[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  num_steps = 1
  dt = ${const_dt}
  petsc_options_iname = '-pc_type -sub_pc_factor_levels -ksp_gmres_restart'
  petsc_options_value = 'lu      2                     200'
  line_search = 'none'
  nl_rel_tol = 1e-4
  nl_abs_tol =  1e-7
  nl_max_its = 60
  automatic_scaling = true
  verbose = true
  off_diagonals_in_auto_scaling = true
[]
[Outputs]
  exodus = true
  print_linear_residuals = false
[]

no AD rather than by hand:

const_num = 1.
const_dt = 1.e-1
# const_scaling = 1e0

vel_x = vel_x
vel_y = vel_y
vel_z = vel_z
pressure = p
[GlobalParams]
  pspg = true
  supg = true
  u = vel_x
  v = vel_y
  w = vel_z
  pressure = p
  integrate_p_by_parts = true
  alpha = 1
[]
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    xmax = ${const_num}
    ymin = 0
    ymax = ${const_num}
    zmin = 0
    zmax = ${const_num}
    nx = 10
    ny = 10
    nz = 10
  []
  [./corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = gen
  [../]
[]

[Variables]
  [./vel_x]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./vel_y]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./vel_z]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./p]
    # scaling = ${const_scaling}
  [../]
[]

[Kernels]
  [mass]
    type = INSMass
    variable = p
  []
  # x-momentum, time
  [x_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_x}
  []
  # x-momentum, space
  [x_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_x}
    component = 0
  []
  # y-momentum, time
  [y_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_y}
  []
  # y-momentum, space
  [y_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_y}
    component = 1
  []
  # z-momentum, time
  [z_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_z}
  []
  # z-momentum, space
  [z_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_z}
    component = 2
  []
[]

[BCs]
  [no_slip_x]
    type = DirichletBC
    variable = ${vel_x}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [no_slip_y]
    type = DirichletBC
    variable = ${vel_y}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [no_slip_z]
    type = DirichletBC
    variable = ${vel_z}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [./lid_x]
    type = DirichletBC
    variable = ${vel_x}
    boundary = 'top'
    value = 1.
  [../]
  [./lid_y]
    type = DirichletBC
    variable = ${vel_y}
    boundary = 'top'
    value = 0.
  [../]
  [./lid_z]
    type = DirichletBC
    variable = ${vel_z}
    boundary = 'top'
    value = 0.
  [../]

  [./pressure_pin]
    type = DirichletBC
    variable = ${pressure}
    boundary = 'pinned_node'
    value = 0
  [../]
[]

[Materials]
  [./const]
    type = GenericConstantMaterial
    prop_names = 'rho      mu               '
    prop_values = '1e4     ${const_num}     '
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  num_steps = 1
  dt = ${const_dt}
  petsc_options_iname = '-pc_type -sub_pc_factor_levels -ksp_gmres_restart'
  petsc_options_value = 'lu       2                     200'
  line_search = 'none'
  nl_rel_tol = 1e-4
  nl_abs_tol =  1e-7
  nl_max_its = 60

  automatic_scaling = true
  verbose = true
  off_diagonals_in_auto_scaling = true
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
[]

# [Debug]
#   show_var_residual_norms = true
# []

The AD and by hand codes in the Navier-stokes module have been checked by me. But I have not found the bugs' reason.

Thanks,

QINGJI SU

[lindsayad]

This hand-coded input version converges about as well as the AD version. The changes were

• scaling grouping for velocity variables (this happens naturally for the vector velocity variable used in AD)
• setting preset = false which applies to the Dirichlet boundary conditions. Vector Dirichlet boundary conditions do not have the preset option, e.g. they are never preset, whereas preset = true is the default for non-vector Dirichlet boundary conditions
• setting transient_term = true in the GlobalParams block. This adds the transient term into the strong residual used for SUPG and PSPG. This happens automatically for AD

const_num = 1.
const_dt = 1.e-1
# const_scaling = 1e0

vel_x = vel_x
vel_y = vel_y
vel_z = vel_z
pressure = p
[GlobalParams]
  pspg = true
  supg = true
  u = vel_x
  v = vel_y
  w = vel_z
  pressure = p
  integrate_p_by_parts = true
  alpha = 1
  transient_term = true
  preset = false
[]
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    xmax = ${const_num}
    ymin = 0
    ymax = ${const_num}
    zmin = 0
    zmax = ${const_num}
    nx = 10
    ny = 10
    nz = 10
  []
  [./corner_node]
    type = ExtraNodesetGenerator
    new_boundary = 'pinned_node'
    nodes = '0'
    input = gen
  [../]
[]

[Variables]
  [./vel_x]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./vel_y]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./vel_z]
    # scaling = ${const_scaling}
    [InitialCondition]
      type = ConstantIC
      value = 1e-15
    []
  [../]
  [./p]
    # scaling = ${const_scaling}
  [../]
[]

[Kernels]
  [mass]
    type = INSMass
    variable = p
  []
  # x-momentum, time
  [x_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_x}
  []
  # x-momentum, space
  [x_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_x}
    component = 0
  []
  # y-momentum, time
  [y_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_y}
  []
  # y-momentum, space
  [y_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_y}
    component = 1
  []
  # z-momentum, time
  [z_momentum_time]
    type = INSMomentumTimeDerivative
    variable = ${vel_z}
  []
  # z-momentum, space
  [z_momentum_space]
    type = INSMomentumLaplaceForm
    variable = ${vel_z}
    component = 2
  []
[]

[BCs]
  [no_slip_x]
    type = DirichletBC
    variable = ${vel_x}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [no_slip_y]
    type = DirichletBC
    variable = ${vel_y}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [no_slip_z]
    type = DirichletBC
    variable = ${vel_z}
    boundary = 'bottom right left front back'
    value = 0.
  []
  [./lid_x]
    type = DirichletBC
    variable = ${vel_x}
    boundary = 'top'
    value = 1.
  [../]
  [./lid_y]
    type = DirichletBC
    variable = ${vel_y}
    boundary = 'top'
    value = 0.
  [../]
  [./lid_z]
    type = DirichletBC
    variable = ${vel_z}
    boundary = 'top'
    value = 0.
  [../]

  [./pressure_pin]
    type = DirichletBC
    variable = ${pressure}
    boundary = 'pinned_node'
    value = 0
  [../]
[]

[Materials]
  [./const]
    type = GenericConstantMaterial
    prop_names = 'rho      mu               '
    prop_values = '1e4     ${const_num}     '
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  num_steps = 1
  dt = ${const_dt}
  petsc_options_iname = '-pc_type -sub_pc_factor_levels -ksp_gmres_restart'
  petsc_options_value = 'lu       2                     200'
  line_search = 'none'
  nl_rel_tol = 1e-4
  nl_abs_tol =  1e-7
  nl_max_its = 60

  automatic_scaling = true
  verbose = true
  off_diagonals_in_auto_scaling = true
  scaling_group_variables = 'vel_x vel_y vel_z; p'
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
[]

# [Debug]
#   show_var_residual_norms = true
# []

[lindsayad]

If I remove the code of automatic_scaling = true, the problem could be solved. That is a strange phenomenon. Therefore I hardly use automatic_scaling = true. But I must admit that in some problems the problem convergence progress becomes better if I use automatic_scaling = true. So I am very tangled if I should use automatic_scaling = true.

Yea automatic scaling is not a silver bullet. Sometimes it is helpful and sometimes it's not. I would suggest always grouping like variables into the scaling computation, like I did in the input above with the scaling_group_variables parameter. If a problem does not solve well with automatic scaling and solve_type = NEWTON, then I would suggest running with -pc_type svd -pc_svd_monitor on a problem with 1000 dofs or less and see whether the problem is well-posed.

[suqingji]

I got it, your answer is helpful to me. Thank you very much.

[lindsayad]

you're welcome! Thanks for your good questions!

[suqingji]

Hi @lindsayad,
Please forgive me for bothering you again.
I encountered a convergence problem when I simulated a flow problem.

If I use the INSAD to simulate:

19 Nonlinear |R| = 8.042014e-16

    |residual|_2 of individual variables:
                   velocity: 6.55322e-16
                   p:        4.66147e-16
  Linear solve did not converge due to DIVERGED_NANORINF iterations 0
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 19

 Solve Did NOT Converge!

  Finished Solving                                                                       [225.63 s] [   81 MB]


Aborting as solve did not converge

Finished Executing

If I use INS to simulate:

60 Nonlinear |R| = 8.810293e-16

    |residual|_2 of individual variables:
                      vel_x: 4.27504e-16
                      vel_y: 4.32017e-16
                      vel_z: 4.25759e-16
                      p:     4.74914e-16
Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations 60

 Solve Did NOT Converge!

  Finished Solving                                                                       [614.75 s] [   84 MB]


Aborting as solve did not converge

Finished Executing

In fact, the output on the screen if I use INS is normal because of my nl_ral_tol , nl_abs_tol and nl_max_its setup. However, the output on the screen, if I use INSAD, is abnormal for Linear solve did not converge due to DIVERGED_NANORINF iterations 0.

As I know, the reason for DIVERGED_NANORINF is that divide by zero when the program solves the Linear equations. That means there are zeros on the diagonal of the matrix. But the same problem if I use the INS rather than INSAD the situations did not happen.

So I suspect that there are some bugs when the AD build the jacobian matrix.

Here are my input files:
INSAD:

[Mesh]
  [./gen]
    type = FileMeshGenerator
    file = dvi-exodus.exo
  [../]
  [inlet_left]
    type = ParsedGenerateSideset
    input = gen
    combinatorial_geometry = 'abs(x + 1500) < 1e-4'
    new_sideset_name = 'inlet_left'
  []
  [inlet_right]
    type = ParsedGenerateSideset
    input = inlet_left
    combinatorial_geometry = 'abs(x - 1500) < 1e-4'
    new_sideset_name = 'inlet_right'
  []
  [./scale]
    type = TransformGenerator
    input = inlet_right
    transform = SCALE
    vector_value ='1e-03 1e-03 1e-03'
  []
[]

[Modules]
  [IncompressibleNavierStokes]
    equation_type = steady-state
    
    velocity_boundary = 'WALL-CORE WALL-DOWN WALL-DVI WALL-HEAT WALL-INLET WALL-NOHEAT WALL-UP inlet_left   inlet_right  INLET-DVI'
    velocity_function = '0 0 0     0 0 0     0 0 0    0 0 0     0 0 0      0 0 0       0 0 0   1 0 0        -1 0 0       4.09 0 -2.87'
    
    initial_velocity = '1e-15 1e-15 1'
    add_standard_velocity_variables_for_ad = false

    density_name = rho
    dynamic_viscosity_name = mu

    use_ad = true
    laplace = true
    family = LAGRANGE
    order = FIRST

    add_temperature_equation = false

    alpha = 1
    supg = true
    pspg = true
  []
[]

[BCs]
  [./outlet]
    type = DirichletBC
    variable = p
    boundary = 'OUTLET'
    value = 0
  [../]
[]
[Debug]
   show_var_residual_norms = true
[]

[Materials]
  [./user_defined]
    type = ADGenericConstantMaterial
    block = FLUID
    prop_names = 'rho mu cp k'
    prop_values = '1 1e-02 4.179e3 0.58 '
  [../]
[]
[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -sub_pc_type -ksp_gmres_restart'
  petsc_options_value = 'bjacobi  asm          200'
  line_search = 'none'
  nl_rel_tol = 1e-20
  nl_max_its = 60
[]
[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]
[Outputs]
  exodus = true
  print_linear_residuals = false
[]

INS:

[Mesh]
  [./gen]
    type = FileMeshGenerator
    file = dvi-exodus.exo
  [../]
  [inlet_left]
    type = ParsedGenerateSideset
    input = gen
    combinatorial_geometry = 'abs(x + 1500) < 1e-4'
    new_sideset_name = 'inlet_left'
  []
  [inlet_right]
    type = ParsedGenerateSideset
    input = inlet_left
    combinatorial_geometry = 'abs(x - 1500) < 1e-4'
    new_sideset_name = 'inlet_right'
  []
  [./scale]
    type = TransformGenerator
    input = inlet_right
    transform = SCALE
    vector_value ='1e-03 1e-03 1e-03'
  []
[]
[Variables]
  [./vel_x]
    [./InitialCondition]
      type = ConstantIC
      value = 1e-15
    [../]
  [../]
  [./vel_y]
    [./InitialCondition]
      type = ConstantIC
      value = 1e-15
    [../]
  [../]
  [./vel_z]
    [./InitialCondition]
      type = ConstantIC
      value = 1
    [../]
  [../]
  [./p]
  [../]
[]
[Kernels]
  # mass
  [./mass]
    type = INSMass
    variable = p
    u = vel_x
    v = vel_y
    w = vel_z
    pressure = p
  [../]
  # x-momentum, space
  [./x_momentum_space]
    type = INSMomentumLaplaceForm
    variable = vel_x
    u = vel_x
    v = vel_y
    w = vel_z
    pressure = p
    component = 0
  [../]
  # y-momentum, space
  [./y_momentum_space]
    type = INSMomentumLaplaceForm
    variable = vel_y
    u = vel_x
    v = vel_y
    w = vel_z
    pressure = p
    component = 1
  [../]
  # z-momentum, space
  [./z_momentum_space]
    type = INSMomentumLaplaceForm
    variable = vel_z
    u = vel_x
    v = vel_y
    w = vel_z
    pressure = p
    component = 2
  [../]
[]

[BCs]
  [./x_no_slip]
    type = DirichletBC
    variable = vel_x
    boundary = 'WALL-CORE WALL-DOWN WALL-DVI WALL-HEAT WALL-INLET WALL-NOHEAT WALL-UP'
    value = 0.0
  [../]
  [./y_no_slip]
    type = DirichletBC
    variable = vel_y
    boundary = 'WALL-CORE WALL-DOWN WALL-DVI WALL-HEAT WALL-INLET WALL-NOHEAT WALL-UP'
    value = 0.0
  [../]
  [./z_no_slip]
    type = DirichletBC
    variable = vel_z
    boundary = 'WALL-CORE WALL-DOWN WALL-DVI WALL-HEAT WALL-INLET WALL-NOHEAT WALL-UP'
    value = 0.0
  [../]

  [./x_in_left]
    type = DirichletBC
    variable = vel_x
    boundary = 'inlet_left'
    value = 1
  [../]
  [./y_in_left]
    type = DirichletBC
    variable = vel_y
    boundary = 'inlet_left'
    value = 0.0
  [../]
  [./z_in_left]
    type = DirichletBC
    variable = vel_z
    boundary = 'inlet_left'
    value = 0.0
  [../]
  
  [./x_in_right]
    type = DirichletBC
    variable = vel_x
    boundary = 'inlet_right'
    value = -1
  [../]
  [./y_in_right]
    type = DirichletBC
    variable = vel_y
    boundary = 'inlet_right'
    value = 0.0
  [../]
  [./z_in_right]
    type = DirichletBC
    variable = vel_z
    boundary = 'inlet_right'
    value = 0.0
  [../]
  
  [./x_in_dvi]
    type = DirichletBC
    variable = vel_x
    boundary = 'INLET-DVI'
    value = 4.09
  [../]
  [./y_in_dvi]
    type = DirichletBC
    variable = vel_y
    boundary = 'INLET-DVI'
    value = 0.0
  [../]
  [./z_in_dvi]
    type = DirichletBC
    variable = vel_z
    boundary = 'INLET-DVI'
    value = -2.87
  [../]
  
  [./outlet]
    type = DirichletBC
    variable = p
    boundary = 'OUTLET'
    value = 0
  [../]
[]
[Debug]
   show_var_residual_norms = true
[]

[Materials]
  [./user_defined]
    type = GenericConstantMaterial
    block = FLUID
    prop_names = 'rho mu cp k'
    prop_values = '1 1e-02 4.179e3 0.58 '
  [../]
[]
[GlobalParams]
  gravity = '0 0 0'
  supg = true
  pspg = true
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -sub_pc_type -ksp_gmres_restart'
  petsc_options_value = 'bjacobi  asm          200'
  line_search = 'none'
  nl_rel_tol = 1e-20
  nl_max_its = 60
[]
[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]
[Outputs]
  exodus = true
  print_linear_residuals = false
[]

Thanks and best wishes.

[suqingji]

Hi, @lindsayad,
First thanks for your reply.
After some tests, I have gotten the reason DIVERGED if I use parallel. The root reason is that the matrix value builds from different PDEs is very different. Therefore, sometimes if we use automatic scaling, the problem maybe can converge. Often we must accord our experience to set a scaling value for those variables. As you said many days ago, automatic scaling is not a silver bullet.

B-U-T, I think there are two methods to solve the problem at the root.

The first one is that change all original PDEs to dimensionless PDEs. So, does your group have the plan to add a dimensionless solver to accelerate the program and improve convergence?

The second one is adding a separate solver which is different from MultiApps. This is what I expect long ago. Of course, I saw a plan that Idaho to add a separate solver to MOOSE at the workshop many months ago. So, could you please introduce some progress about the work?

Best wishes.

[lindsayad]

So, does your group have the plan to add a dimensionless solver to accelerate the program and improve convergence?

No we don't. Note that the MOOSE framework, libMesh, and PETSc have no notion of dimensions. It's up to the applications (or modules) how they want to handle units. The primary developers of the navier_stokes module, @GiudGiud and me, have no funding mandate to work on units. We always welcome contributions from the community, however.

The second one is adding a separate solver which is different from MultiApps. This is what I expect long ago. Of course, I saw a plan that Idaho to add a separate solver to MOOSE at the workshop many months ago. So, could you please introduce some progress about the work?

We merged #22226 not very long ago that allows multiple nonlinear systems in the same input file. You can create your own custom Executioner to leverage this capability if you wish. The MOOSE team is also developing some TH specific segregated solve schemes outlined in #22356 and (partially tackled) in #22699

[suqingji]

Ok, thanks for your reply.

Sorry, Maybe my description makes you misunderstand my means. A dimensionless solver(https://en.wikipedia.org/wiki/Dimensionless_quantity) is not a solver which is irrelevant to dimension(https://en.wikipedia.org/wiki/Dimension).

If your group does have not a plan to add the dimensionless solver to MOOSE, I will add a plan to my plan list to contribute the one to MOOSE.

Secondly, I appreciate your group for https://github.com/idaholab/moose/pull/22226.

Best wishes.

[lindsayad]

I'm well acquainted with dimensionless quantities. From MOOSE's perspective everything is just a number so from that perspective we are always dimensionless. It is the user and the application which generally applies dimensional meanings in their inputs, material properties, etc. For instance an application could choose to code up the Navier-Stokes equations using just the Reynolds number ... but the navier_stokes developers did not make that choice

[lindsayad]

If your group does have not a plan to add the dimensionless solver to MOOSE, I will add a plan to my plan list to contribute the one to MOOSE.

We always love contributions!

[suqingji]

Hi @lindsayad
Please forgive me for bothering you again.

When I debugged I found a problem in MOOSE. Now please admit me introduce the problem to you.

Firstly, let us consider a transient heat conduction problem. Here is the PDE of the heat conduction problem.

If we consider the problem as one dimension. Set the params as:
kappa = 0.59
rho = 1000
cp = 4.2e3
Lx = 1
dx = 0.01
nx = 100
dt = 10
And Set boundary conditions as:
T_left = 1
T_right = 0

Now, let us create an input file. Here are the details of the input file.

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 100
  []
[]

[Variables]
  [./T]
  [../]
[]

[Kernels]
  [./heat]
    type = ADHeatConduction
    variable = T
  [../]
  [./heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = T
  [../]
[]

[Materials]
  [./const]
    type = ADGenericConstantMaterial
    prop_names = 'density   specific_heat     thermal_conductivity'
    prop_values = '1000     4.2e3             0.59'
  [../]
[]

[BCs]
 [./left]
   type = DirichletBC
   variable = T
   boundary = left
   value = 1
 [../]
 [./right]
   type = DirichletBC
   variable = T
   boundary = right
   value = 0
 [../]
[]
[Preconditioning]
  [./fdp]
    type = FDP
    full = true
  [../]
[]
[Executioner]
  type = Transient
  solve_type = 'NEWTON'

  dt = 10
  end_time = 10
[]

[Outputs]
  exodus = true
[]

The solution is shown here:

Obviously, there are some values that are bigger than 1 and smaller than 0.

In fact, this is an abnormal phenomenon.

Secondly, let me prove my comment.
I write a program for the question through the libMesh which also is the MOOSE's underlying framework.

This is the plot calculated by libMesh:

At last, I use python to write a program for the question through the finite difference method
This is the plot calculated by FDM with Python:

So, I suspect that there are some bugs in the Time Derivative part in MOOSE.

Thanks.

[suqingji]

Hi, @lindsayad, @GiudGiud, @YaqiWang

As I mentioned above, in the transient diffusion problem there is an oscillation problem. This is caused by CG FEM's non-local conservation.

Fortunately, there is a method that can be used to solve the problem. The method is lumping the transient term. I also found that in the PorousFlow module there is a class that has implemented the lumped transient term. Here is the detail page.

But unfortunately, in the Navier-Stokes module, if I want to couple a variable using some codes such as &adCoupledDot("temperature"), there is no lumped coupled method. I try to add the method to Base Class, but I failed.

Would you like to add the method to the class Coupleable.

Thanks. Best wishes.

[GiudGiud]

Hello

Which method are you missing in Coupleable? adCoupledDot does exist
The massLumpedTimeDerivative would likely only be added to NS if useful, coupleable is a very general framework-wide class

Guillaume

[suqingji]

Hi,

I want to add a "lumped" transient term of the Momentum equation and Energy Equation.

As you know, in the INSAD module, all the terms are computed in the Material Class. So if I want to compute a lumped transient term, I must couple the variable as a lumping method in the Material Class. But I did not find any lumping coupled variable method, especially vector variable.

Thanks.

[suqingji]

Hello,

For example, in non-AD type:

in AD type:

So, I think this is a needed method.

Thanks, and best wishes.

[lindsayad]

Yes we would be interested in having that implemented. You can create an issue for this, and I can probably get to it in the not too distant future. It shouldn't take much time.