Mortar-Based Mechanical Contact Performance

Results of performance studies on the mortar-based contact approach are shown here.

Frictionless contact algorithm comparison

ConstraintDisplacementNCP functionTime (arbitrary units)Time stepsNonlinear iterations
NodalMortarMin4.16440104
NodalMortarFB5.02040135
NodalNodalMin3.12441104
NodalNodalFB4.01441149
MortarMortarMin4.46140106
MortarMortarFB5.57740136
NodalNodalRANFS2.7004099

The first column denotes the discretization algorithm used for applying the frictionless contact constraints. Nodal denotes use of a NodeFaceConstraint; Mortar denotes use of a MortarConstraint. The second column denotes the discretization used for applying the contact forces to the displacement residuals. The third column denotes the type of non-linear complementarity problem (NCP) function used to ensure that the contact constraints are satisfied. Min indicates the canonical min function (see std::min); FB represents the Fischer-Burmeister function. RANFS denotes the Reduced Active Nonlinear Function Set scheme in which no Lagrange Multipliers are used, and instead the non-linear residual equations at the secondary nodes are replaced with the gap function. The fourth column in the table is the simulation time in arbitrary units (since timings will be different across machines). The fifth column is the number of time steps required to reach the simulation end time. The final, sixth column is the cumulative number of non-linear iterations taken during the simulation (note that this does not include any non-linear iterations from failed time steps).

Notes:

  • Clearly having mortar mesh generation slows the simulation down, which is not surprising

  • The min NCP function is undeniably better for solving normal contact

  • For the pure nodal algorithms, the time step that did not converge featured classic ping-ponging behavior:


 5 Nonlinear |R| = 4.007951e-04
    |residual|_2 of individual variables:
                  disp_x:    0.000399808
                  disp_y:    2.75599e-05
                  normal_lm: 5.52166e-06


The number of nodes in contact is 11

      0 Linear |R| = 4.007951e-04
      1 Linear |R| = 1.287307e-04
      2 Linear |R| = 8.423398e-06
      3 Linear |R| = 1.046825e-07
      4 Linear |R| = 8.017310e-09
      5 Linear |R| = 3.053040e-10
  Linear solve converged due to CONVERGED_RTOL iterations 5
 6 Nonlinear |R| = 4.432193e-04
    |residual|_2 of individual variables:
                  disp_x:    0.000396694
                  disp_y:    0.00019545
                  normal_lm: 2.96013e-05


The number of nodes in contact is 11

      0 Linear |R| = 4.432193e-04
      1 Linear |R| = 1.355935e-04
      2 Linear |R| = 1.216010e-05
      3 Linear |R| = 6.386952e-07
      4 Linear |R| = 2.235594e-08
      5 Linear |R| = 2.884193e-10
  Linear solve converged due to CONVERGED_RTOL iterations 5
 7 Nonlinear |R| = 4.008045e-04
    |residual|_2 of individual variables:
                  disp_x:    0.000399816
                  disp_y:    2.76329e-05
                  normal_lm: 5.29313e-06


The number of nodes in contact is 11

      0 Linear |R| = 4.008045e-04
      1 Linear |R| = 1.287272e-04
      2 Linear |R| = 8.423081e-06
      3 Linear |R| = 1.047782e-07
      4 Linear |R| = 8.054781e-09
      5 Linear |R| = 3.046073e-10
  Linear solve converged due to CONVERGED_RTOL iterations 5
 8 Nonlinear |R| = 4.432194e-04

Frictional contact algorithm comparison

LM normalLM tangentialDisplacementNCP function normalNCP function tangentialTime (arbitrary units)Time stepsNonlinear iterationsCLI PETSc options
MortarMortarMortarFBFB8.24140175None
MortarMortarMortarMinFB7.92840159None
NodalMortarMortarMinFB7.45940152None
MortarMortarMortarMinMin11.23741234None
NodalNodalMortarMinMin39.40955275-snes_ksp_ew 0
NodalNodalMortarFBFBNANANANone

Notes:

  • NA: solve did not converge

  • Timings run on a different machine than the frictionless cases

  • The most performant case uses a NodeFaceConstraint discretization for enforcing the normal contact conditions and MortarConstraint discretizations for enforcement of the Coulomb frictional constraints and application of forces to the displacement residuals. Interestingly, this performant case uses different NCP functions for normal and tangential constraints: std::min for the former and Fischer-Burmeister for the latter. This performant case is used for comparison with the node-face penalty algorithm, shown below:

NCP-LM-Mortar vs Penalty-NodeFace

The table below compares the timing and solver performance of NCP-LM-Mortar and Penalty-NodeFace algorithms. NCP-LM refers to use of an NCP function for contact constraint enforcement on a lagrange multiplier. The "Mortar" designation denotes that a mortar discretization is used for enforcing the tangential Coulomb friction conditions and applying contact forces to the displacement residuals.

AlgorithmTime (arb. units)Time steps to end timeCumulative non-linear iterations
NCP-LM-Mortar13.901151476
Penalty-NodeFace20.711151938

There's a cost associated with generation of the mortar segment mesh that partially offsets the fact that the mortar case takes nearly half the non-linear iterations of the penalty case.

Petsc options for contact

Recommended PETSc options for use with mortar based frictional contact are:


[Executioner]
  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
  petsc_options_iname = '-pc_type -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu       1e-5          NONZERO               1e-15'
[]

Using Eisenstat-Walker is advantageous for frictional contact because time is not wasted in the linear solve in early non-linear iterations while the contact set and stick/slip conditions are being resolved. Later in the non-linear solve when the set of constraints has been resolved, more linear iterations will be used as the non-linear solver moves through the quadratic basin. Experience has shown that a choice of 1e-5 for the matrix free finite differencing parameter works well for many problems. However, the user may want to experiment with values anywhere between 1e-8 and 1e-4 depending on their multi-physics. A very small non-zero shift is used to avoid zero pivots during the LU decomposition. This may be extraneous in many cases. Note that the Jacobian entries for mortar based contact are accurate and complete enough that incomplete factorization may be used in serial or as a sub-block solver for block jacobi or additive schwarz in parallel. This may be necessary for large problems where lu does not scale.

The recommended PETSc options for use with NodeFaceConstraint based contact are shown below :


[Executioner]
  ...
  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap
                        -ksp_gmres_restart'
  petsc_options_value = 'asm lu 20 101'
  ...
[../]

Scaling effects

The effects of scaling on non-linear solve convergence are shown below for different kinds of contact formulations. The results are based on the input files for RANFS, kinematic, and tangential penalty, running with two processes.

SchemeCumulative nonlinear iterationsCumulative linear iterationsInitial condition numberVariable Scaling Factor
SMP PJFNK RANFS no scaling AMG685109e31
SMP PJFNK RANFS residual auto-scaling AMG653911e21.1e-2
SMP PJFNK RANFS Jacobian auto-scaling AMG663724e14.1e-4
SMP PJFNK Kinematic no scaling AMG583059e31
SMP PJFNK Kinematic residual auto-scaling AMG583051e21.1e-2
SMP PJFNK Kinematic Jacobian auto-scaling AMG583054e14.1e-4
SMP PJFNK Tangential Penalty no scaling AMG654009e31
SMP PJFNK Tangential Penalty residual auto-scaling AMG654001e21.1e-2
SMP PJFNK Tangential Penalty Jacobian auto-scaling AMG654004e14.1e-4
FD PJFNK RANFS jacobian auto-scaling LU54704e14.1e-4
FD PJFNK Kinematic jacobian auto-scaling LU62954e14.1e-4

Important takeaways:

  • The solve efficiency of kinematic and tangential penalty formulations is independent of scaling (within the window of this problem)

  • Because RANFS constraint residuals/Jacobians are of gap magnitude, the RANFS solve really does perform better when the internal physics is scaled to be on the same order of magnitude

  • There are some bugs with the RANFS Jacobian functions because RANFS takes more nonlinear and linear iterations than kinematic, but with a FD Jacobian RANFS takes significantly less nonlinear and linear iterations than kinematic