Function Jacobian Pairs

[JaredCrean2]

Overview

Currently the code passes around functions (for example, Newton's method takes the function it is solving as an argument). This also allows computing the Jacobian via coloring. When computing the Jacobian explicitly (not via coloring), the Jacobian of the physics is calculated. This raises the possibility that explicit Jacobian calculation could produce a different Jacobian than coloring if Newton's method is passed a non-physics function.

The only occurs if the function physicsJac is used to compute the Jacobian of something other than a physics. Arguably this is an incorrect usage of physicsJac, however it is rather useful for methods like homotopy, which constructs a fictitious physics.

Currently this bug only affects the homotopy, and it can be easily worked around. Eventually a more through solution should be developed.

Technical details

If the ctx_residual[1] passed into newtonInner is not a physics function and this ctx_residual gets passed to physicsJac, the bug will occur.

Possible Solutions

I think think of two solutions:

1. Instead of passing around functions, pass around an object that contains both both the function and its Jacobian. ie.

   type FunctionJacobianPair
     func::Function  # evaluates the function
     jac::Function  # evaluates the Jacobian of the function
   end

2. Require any method that wants to create a fictitious physics to create an AbstractSolutionData to go with it and extend the proper methods (including the Jacobian calculation function) for it.

[JaredCrean2]

This problem came up again when trying to write a block Jacobi smoother. To make the smoother general, I need to do matrix-free Jacobian vector products with nonlinear problem Jacobian.