Variational inequality solver and non-linear elements

[TheGreatCid]

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

Hi all,

There appears to be an issue with petsc's inequality solver and non-linear elements: The vinewtonrsls petsc method seems to be incorrectly enforcing constraints.

For example, the following input file has a ADTimeDerivative and ADBodyForce kernel. The body force kernel is trying to drive the variable below zero, however, there is a constraint to not let the variable go below zero.

Attached are the mesh files for the tet10's and tet4's
box_tet.zip

[Mesh]
    [fmg]
        type = FileMeshGenerator
        file = 'box_tet10.msh'
    []
[]

[Variables]
    [d]
        order = SECOND
    []
[]

[AuxVariables]
    [bounds_dummy]
        order = SECOND
    []
[]

[Bounds]
    [irreversibility]
        type = VariableOldValueBounds
        variable = bounds_dummy
        bounded_variable = d
        bound_type = lower
    []
[]

[Kernels]
    [diff]
        type = ADTimeDerivative
        variable = d
    []
    [body]
        type = ADBodyForce
        variable = d
        value = -0.1
    []
[]

[Executioner]
    type = Transient
    dt = 1
    end_time = 2
    solve_type = 'NEWTON'

    petsc_options_iname = '-pc_type -snes_type -snes_vi_monitor'
    petsc_options_value = 'lu vinewtonrsls true'
    nl_abs_tol = 1e-08
    nl_rel_tol = 1e-06


     [Quadrature]
         order = FOURTH
     []
[]
[Outputs]
    exodus = true
[]

When using linear elements, this is enforced correctly (d=0 on all nodes), however, when using non-linear tet10's, the following field is generated:

I don't have any leads for this. I checked the bounds vector being generated by BoundsBase.C and it seems correct.

Any insight would be appreciated!

[GiudGiud]

hi

thanks for reporting. I ll take a look

[GiudGiud]

I can confirm the problem.
We need to figure out why the constraints are being reduced from 63 constraints to 49 (which seems to be the issue)

[maxnezdyur]

Applying the force in a kernel did work.
Works with hex20, tet10, tet14 with second order.
Results:
Max absolute value of the d variable.
The VCP works well with it too.

Postprocessor Values:
+----------------+----------------+
| time           | max_d          |
+----------------+----------------+
|   0.000000e+00 |   0.000000e+00 |
|   1.000000e+00 |   3.388132e-20 |
|   2.000000e+00 |   3.388132e-20 |
|   3.000000e+00 |   3.388132e-20 |
+----------------+----------------+

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1.0
    zmin = 0
    zmax = 1
    nx = 2
    ny = 2
    nz = 2
    elem_type = TET10
  []
[]

[Variables]
  [d]
    order = SECOND
  []
  [lm]
    order = SECOND
  []
[]

[NodalKernels]
  [positive_constraint]
    type = LowerBoundNodalKernel
    variable = lm
    v = d
    lower_bound = 0
  []

[]

[Kernels]
  [diff]
    type = ADTimeDerivative
    variable = d
  []
  [body]
    type = ADBodyForce
    variable = d
    function = '-1 * abs(sin(x)+cos(y) + sin(y)*cos(x))'
  []
  [forces]
    type = CoupledForce
    variable = d
    v = lm
  []

[]
[Preconditioning]
  [vcp]
    type = VCP
    lm_variable = lm
    primary_variable = d
    preconditioner = 'AMG'
  []
[]

[Executioner]
  type = Transient
  dt = 1
  end_time = 3
  solve_type = 'NEWTON'

  # petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu                        NONZERO
  # 1e-15'

  nl_abs_tol = 1e-08
  nl_rel_tol = 1e-06
  num_steps = 3
[]

[Postprocessors]
  [max_d]
    type = NodalExtremeValue
    variable = d
    value_type = max_abs
  []

[]
[Outputs]
  exodus = true
[]

[lindsayad]

That's great! Still would be best to get the bounds deviation weighted, but applying the force in a variational way is definitely a good improvement, and good to hear it yielded results

[maxnezdyur]

This works, is this what you meant?

Kernels]
  [diff]
    type = ADTimeDerivative
    variable = d
  []
  [body]
    type = ADBodyForce
    variable = d
    function = '-1 * abs(sin(x)+cos(y) + sin(y)*cos(x))'
  []
  [forces]
    type = CoupledForce
    variable = d
    v = lm
  []
  [parsed_lowerbound]
    type = ADMatBodyForce
    variable = lm
    material_property = lower_bound
  []

[]
[Materials]
  [parsed_lower_bound]
    type = ADParsedMaterial
    property_name = 'lower_bound'
    expression = 'min(lm,d-0) #change 0 to any constant'
    coupled_variables = 'd lm'
  []
[]

Final input file: Implements irreversability constraint

[Mesh]
  [MeshGenerator]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1.0
    zmin = 0
    zmax = 1
    nx = 2
    ny = 2
    nz = 2
    elem_type = TET14
  []
[]

[Variables]
  [d]
    order = SECOND
    initial_condition = 1.0
  []
  [lm]
    order = SECOND
  []
[]

[Kernels]
  [diff]
    type = ADTimeDerivative
    variable = d
  []
  [body]
    type = ADBodyForce
    variable = d
    function = '-1 * abs(sin(x)+cos(y) + sin(y)*cos(x))'
  []
  [forces]
    type = CoupledForce
    variable = d
    v = lm
  []
  [parsed_lowerbound]
    type = ADMatBodyForce
    variable = lm
    material_property = lower_bound
  []

[]
[Materials]
  [parsed_lower_bound]
    type = ADDerivativeParsedMaterial
    property_name = 'lower_bound'
    expression = 'min(lm,d-d_old)'
    coupled_variables = ' d lm'
    material_property_names = 'd_old:=Old[d_mat]'
  []
  [d_mat]
    type = ParsedMaterial
    property_name = 'd_mat'
    coupled_variables = d
    expression = 'd'
  []
[]
[Preconditioning]
  [vcp]
    type = VCP
    lm_variable = lm
    primary_variable = d
    preconditioner = 'AMG'
  []
[]

[Executioner]
  type = Transient
  dt = 1
  end_time = 3
  solve_type = 'NEWTON'

  # petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
  # petsc_options_value = 'lu                        NONZERO
  # 1e-15'

  nl_abs_tol = 1e-08
  nl_rel_tol = 1e-06
  num_steps = 3
[]

[Postprocessors]
  [max_d]
    type = NodalExtremeValue
    variable = d
    value_type = max_abs
  []

[]
[Outputs]
  exodus = true
[]

[lindsayad]

It's close. But you don't integrate the min function. You're supposed to integrate d - foo against the test function alone, and then perform the min comparison with the integrated constraint vs. the associated LM dof value

[TheGreatCid]

@lindsayad Hi! I just want to revive this issue to see if there has been any discussion on something similar elsewhere. That is getting a the petsc variational solver to work with TET10 elements.