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

[Smith, Barry F.]

  Alex,

     This is the first step in understanding the situation and far from the last.

      I would start by trying to think about the process more abstractly without thinking about post-solves, pre-solves etc since those details and their mechanics often get in the way of coming up with a good solution. In other words don't even think about it in terms of PETSc.

      So I start by thinking about 

            Outer loop - the contact problem 
                  Inner loop - try to solve the nonlinear problem with the given contacts
                  Do Newton until 
                         either satisfactory reduction in residual or
                         residual is stuck and one needs to change the contact. 
                   check if contacts should be changed (note this has two flavors 
                                                                              (1) the Newton got stuck or 
                                                                              (2) Newton did not get stuck but one can improve the residual by changing a contact

          I would never update the contact state inside Newton, inside a function evaluation or inside a line search, it should be done outside of all of those things. This is how one solves DVI problems.

   Let's iterate now by email until we understand each other,

   Barry




> On Dec 20, 2017, at 8:59 AM, Alexander Lindsay <alexlindsay239 at gmail.com> wrote:
> 
> 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
>