FEAP interface to PETSc SNES

[Colin McAuliffe]

Hello,

I am trying to make a FEAP interface to SNES in a user macro. The SNES object requires a
function to evaluate the jacobian, where the jacobian Mat is passed as a dummy argument.

It would be nice to be able to call formfe to assemble the jacobian, but the convention
with formfe is to store the jacobain Mat in a common variable. Is there a way to resolve
this conflict without having to write a version of formfe where thejacobian is passed as
a dummy variable? Also, I would like to avoid copying the jacobian to a new Mat if possible.

Thanks
Colin

[FEAP_Admin]

What conflict are you referring to?  For SNES I thought you just need two cover functions, one to set the residual and one for the tangent.  You should be able to construct two cover functions that call formfe.f for you to do this.  The fact that the pointer to the tangent is in a common block should not matter.  Your cover function can be set up to mass the pointer value out as one of its arguments.

[Colin McAuliffe]

Ok I see what you mean. Originally I had Kmat as a dummy which conflicted with the Kmat in common. This seems to be working now with the Kmat pointer being copied and passed through.

Thanks
Colin

[FEAP_Admin]

Great that you have it working.  If you have a brief outline of how you have set up SNES, it would great if you could make a posting to guide others.

[Colin McAuliffe]

Well after some more testing I spoke too soon saying that it is working. Using newton with no line search as test, I am able to converge to meaningful answers with snes but convergence is much slower than with feap's newton solver so I must have done something wrong at some point. When I figure this out I'll post a summary of the steps, and I'd also be willing to post the umacr, if that is permitted on these forums.

[Colin McAuliffe]

Here is a basic outline of what is needed to make a FEAP umacr that PETSc SNES:

First, PETSc SNES requires 2 basic user-provided routines, one routine to evaluate the residual and one routine to assemble the jacobian. Say these are named FormFunction and FormJacobian. Both take the current iterate computed by SNES as an argument. Adding these routines to the SNES context is done with the functions SNESSetFunction and SNESSetJacobian. Now, when SNES needs to evaluate the residual or jacobian, it will call FormFunction and FormJacobian.

FormFunction will go something like this:

1. compute the increment between SNES current iterate x (passed as an argument by SNES) and FEAP current solution u.
2. save FEAP solution and rate arrays in temporary storage
3. update FEAP solution and rate arrays using the increment computed in 1. Note, history variable should not be updated at this point
4. call formfe with isw = 6
5. copy the feap residual array into the PETSc Vec which is to be the output. Note, the FEAP and SNES residuals have opposite signs.
6. restore feap solution and rate arrays from temporary storage

FormJacobian will essentially be the same but with isw = 3 in step 4. Also, if the PETSc ON command has been used formfe will assemble the jacobian to a PETSc Mat called Kmat. The value of this pointer can be passed out of FormJacobian. If a preconditioner is used, it will be computed in this routine as well.

The last ingredient is a way to update FEAP's solution, rate, and history variables after the SNES nonlinear method accepts a solution step. This can be done by creating a custom monitor, which should do steps 1,3,4 above, but with the flag for history variable update turned on. The monitor can also be used to write information to the feap log and output files, and/or to the screen.

The problem I mentioned in the previous post was due to my FormFunction and FormJacobian routines calling SNESGetSolutionUpdate and using that to do the update in step 3. This is not equivalent to the increment computed by the method in step 1.

I hope this is helpful!

Colin

[FEAP_Admin]

Thanks Colin. 
Can you say anything yet about the performance of SNES versus FEAP's built-in Newton/Quasi-Newton methods?  Are there important advantageous to using SNES?

(fyi, it is fine to post your own user code here as well as snippets from FEAP itself)

[Colin McAuliffe]

The advantage of using SNES over FEAP would be if you have a problem where using special nonlinear solver like Jacobian Free Newton Krylov (JFNK) would be useful. For example, if you have a very large problem, JFNK can save memory since you do not need to store the assembled Jacobian matrix; only the preconditioner is assembled. This can also be used for user elements which have not been programmed to compute the jacobian since only residual evaluations are neccesary. Full Approximation Storage (FAS) nonlinear multigrid solvers are also available.

On the other hand, using Newton with no line search or trust region method (-snes_ls basic) will just give you a slower version of FEAP's Newton, since the SNES FormFunction and FormJacobian require some redundant calculations. I attached some log files of example 8 from the example manual using FEAP's Newton and SNES Newton with default tolerances (Relative Energy of step < 1e-16 for feap and relative residual < 1e-8 for SNES) and FEAP's Newton is still faster despite the stricter tolerance requiring more linear solves. I don't know if SNES cubic line search or trust region offer any advantages over FEAP's line search.

Note, the umacr I attached is the minimum that is needed to access SNES from FEAP. For JFNK to really be useful, the routine will have to be modified to compute a preconditioner, and to tell PETSc the sparsity structure of the Jacobian. Also, I have only tested this on one processor.

To use the umacr you would replace

PETSc,ON
LOOP,,NITER
  TANG,,1
NEXT

with

PETSc,ON
SNES