[petsc-users] Mechanical contact -> Changing form of residual function

[Alexander Lindsay]

This question comes from modeling mechanical contact with MOOSE; from
talking with Derek Gaston, this has been a topic of conversation before...

With contact, our residual function is not continuous. Depending on the
values of our displacements, we may or may not have mechanical contact
resulting in extra residual terms from contact forces. We currently have
the following non-linear algorithm, which is just wrong:

Sometime past the first initial non-linear residual evaluation...

Update contact state to state n
Evaluate Jacobian (contact state n)
KSPSolve with Jacobian in contact state n, residual in contact state (n - 1)
Apply line search; evaluate residual in contact state n

Back to top of non-linear loop

Update contact state to state (n + 1)
Evaluate Jacobian (contact state n + 1)
KSPSolve with Jacobian in contact state (n+1), residual in contact state (n)
Apply line search; evaluate residual in contact state (n+1)

...

So clearly this is stupid because we're conducting a linear solve with
J(-du) = R, with R = R_n(u) and J = dR_{n+1}(u) / du, e.g. the Jacobian is
the derivative of our function in an updated contact state whereas our
actual function is in an older state.

What's the right thing to do here? We could update the contact state before
evaluating the residual, e.g. make the contact update part of our residual
evaluation routine called by SNESComputeFunction. However, we often run
using MFFD, which means there's a chance that when evaluating our Jacobian
action within the linear solve, we might cross the discontinuity when
perturbing the residual.

In my mind, a good option would be to update the contact state only when
evaluating non-linear residuals, e.g. when evaluating the initial
non-linear residual or when applying the line search. Is there a hook into
Petsc for doing such a thing? Alternatively, I see that there's an update
function called at the beginning of Petsc's non-linear iteration that can
be set with SNESSetUpdate. There we could update our contact state and then
re-evaluate F using the updated contact state. This would have the effect
that our last evaluated non-linear residual which was tested for
convergence would be different from the new residual that we are inserting
into the RHS of our new KSPSolve, e.g. the functions were evaluated at
different contact states using the same non-linear solution vector.

Anyways, this is a long message and if it's not initially clear, I
apologize in advance. We're aware that Newton can struggle with
discontinuous functions so we're already in a tough spot, but we'd like to
be as algorithmically correct as possible within the Newton context.

Alex
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171220/43d0b37f/attachment.html>