Jacobian calculation

Special mention of calcJacobianComplex.

One of the most involved parts of Newton's method is forming the Jacobian. The NonlinearSolvers module contains the function physicsJac to compute the Jacobian of the physics $\frac{\partial R(q)}{\partial q}$, where $R$ is evalResidual. This function can be used by other methods as a starting point for computing the residual of $g(R(q))$ described in Newton's method page.

NonlinearSolvers.physicsJac — Function.

This function computes the Jacobian of an evalResidual-like function. Specifically, it computes \partial (eqn.q_vec)/ \partial (eqn.res_vec).

Other users of newtonInner (for example, implicit time marching methods) will need to implemnt their own version of this function. The documentation for this function describes the requirements for the other implementations.

`newtonInner guarantees that eqn.q and eqn.q_vec will be consistent when this function is called. When using finite differences, eqn.res and eqn.res_vec must contain the residual of the physics.

Inputs:

mesh: an AbstractMesh
sbp: an SBP operator
eqn: an AbstractSolutionData, eqn.res and eqn.res_vec may be overwritten
opts: options dictonary
jac: the Jacobian, can be an Array, SparseMatrixCSC, or PetscMat
ctx_residual: a tuple of values. ctx_residual[1] must be an evalResidual-like function (eqn.q -> eqn.res) function with signature func(mesh, sbp, eqn, opts, t). See physicsRhs for a more thorough description of ctx_residual.
t: simulation time

Options Keys:

jac_type
jac_method
epsilon
calc_jac_explicit

Implementation Notes:

This function should not have to do any parallel communication. newtonInner ensures that the rhs_func is called before jac_func, and rhs_func handles the parallel communication.

Implementations of this function may perform either (eqn.q -> jacobian) or (eqn.q_vec -> jacobian). The first may be more computationally efficient, but the second can be simpler for debugging.

This function supportes several types of jacobians (dense arrays, SparseMatrixCSC, PetscMat), and several methods for calculating them (finite difference and complex step). Any function calling this function should support them as well.

It is strongly recommneded to use this function to compute the spatial jacobian and them modify the resulting matrix (this function zeros the Jacobian matrix)

When using Petsc matrices, the function may do intermediate assemblies (PETSC_FLUSH_ASSEMBLY), but does not need to do the final assembly.

physicsJac calls one of several functions to compute the Jacobian. Users should not call these functions directly, they should use physicsJac. There are two approaches to computing the Jacobian: using finite differences/complex step and calculating the entries of the matrix explicitly. The latter approach is much faster than the former, however the former does not require the physics module differentiate evalResidual.

Finite Differencing/Complex Step

The functions calcJacFD and calcJacobianComplex compute the Jacobian as a dense matrix, one column at a time. This is very slow and only useful for debugging small cases. A more efficient method is implemented in calcJacobianSparse, which uses a graph coloring approach to perturb the solution at several nodes simultaneously.

NonlinearSolvers.calcJacFD — Function.

NonlinearSolvers.calcJacFD

This function calculates the Jacobian using finite differences, perturbing one degree of freedom at a time. This is slow and not very accurate. The Jacobian is calculated about the point in eqn.q_vec.

Inputs: mesh: AbstractMesh sbp: SBP operator eqn: AbstractEquation object opts: options dictionary pert: perturbation to use the finite differences. Must be of type Tsol. func: residual evaluation function res_0: vector containing residual at the point the Jacobian is calculated

Inputs/Outputs: jac:: Jacobian matrix to be populated. Must be a dense matrix

Aliasing restrictions: res_0 must not alias eqn.res_vec

At the start, calcJacFD assumes: The Jacobian will be calculated at the state that is specified in eqn.q_vec . res_0 should have the residual at that state in it

At exit, eqn.q_vec will have the same values as at the start.

eqn.q and eqn.res will be overwritten in the course of this function.

NonlinearSolvers.calcJacobianComplex — Function.

NonlinearSolvers.calcJacComplex

This function calculates the Jacobian (dense) using the complex step method, perturbing one degree of freedom at a time. This is very slow. The jacobian is calculated about the point in eqn.q_vec.

Inputs: mesh: AbstractMesh sbp: SBP operator eqn: AbstractEquation object opts: options dictionary pert: perturbation to use. Must be of type Tsol. func: residual evaluation function

Inputs/Outputs: jac:: Jacobian matrix to be populated. Must be a dense matrix

Aliasing restrictions: res_0 must not alias eqn.res_vec

NonlinearSolvers.calcJacobianSparse — Function.

NonlinearSolvers.calcJacobianSparse

This function calculate the Jacobian sparsely (only the entries within the sparsity bounds), using either finite differences or algorithmic differentiation. The jacobian is calculated about the point stored in eqn.q (not eqn.q_vec). A mesh coloring approach is used to compute the jacobian. Both eqn.q and the MPI send and receive buffers are perturbed during this process.

Inputs: mesh: AbstractMesh sbp: SBP operator eqn: AbstractEquation object opts: options dictionary pert: perturbation to use for the algorithmic differentiation. Currently, only complex numbers are supported. func: residual evaluation function (eqn.q -> eqn.res) res_0: element-based (3 dimensional) array containing the residual evaluated at the point where the Jacobian is being calculated. This is only used for finite differences (can be a 0 x 0 x 0 array otherwise).

Inputs/Outputs: jac: Jacobian matrix. Must be a sparse matrix type of some kind, (including PetscMat).

Aliasing restrictions: res_0 must not alias eqn.res

NonlinearSolvers.applyPerturbation — Function.

NonlinearSolvers.applyPerturbation

This function applies a perturbation to a the specified degree of freedom on each element according to a mask.

Because this is element based perturbation, opts["parallel_data"] must be "element".

Inputs: mesh: an AbstractMesh color: the color to perturb pert: perturbation to apply. Can be any datatype i: local degree of freedom number (in range 1:numDofPerNode) to perturb j: local node number (in range 1:numNodesPerElement) to perturb

Inputs/Outputs: arr: element based (3D) array of values to perturb shared_data: array of SharedFaceData for ghost elements to be perturbed The receive buffers are perturbed according to the masks in mesh.shared_element_colormasks, the send buffers are perturbed consistently with arr.

Aliasing restrictions: none

NonlinearSolvers.assembleElement — Function.

Assembles the volume integral contribution.

This function assembles the jacobian contribution for a single element into the matrix, when the jacobian is computed using finite differences. This is used by coloring-based methods for computing the jacobian.

Inputs:

helper: a AssembleData object
mesh: AbstractMesh object
eqn: AbstractEquation
res_arr: element-based (3D) array of perturbed residual values
res_0: element-based (3D) array of non-perturbed residual values
el_res: element number of the element we are observing the change in
el_pert: element number of the element that was perturbed
dof_pert: the degree of freedom number of the perturbed dof
epsilon: magnitude of perturbation

Inputs/Outputs:

jac: any kind of matrix (dense, SparseMatrixCSC, PetscMat)

Aliasing restrictions: res_arr and res_0 must not alias each other.

Same as other method, but for complex numbers. See that method for details. res_0 is not used in this case

Inputs

helper: an _AssembleElementData
mesh: a mesh
elnum: element number
jac: 4 dimensional array containing the jacobian of the element

jac contains the data for the jacobian of the volume terms for a given element. $jac[i, j, p, q] = \partial R[i, p, elnum] / \partial eqn.q[j, q, elnum]$. Its size is numDofPerNode x numDofPerNode x numNodesPerElement x numNodesPerElement.

NonlinearSolvers.calcJacCol — Function.

NonlinearSolvers.calcJacCol

This function extracts the entries for one column of the Jacobian from two residual evaluates that come from finite differences.

Inputs: res_0: vector of unperturbed residual values res: vector of perturbed residual values epsilon: magnitude of perturbation

Inputs/Outputs: jac_row = vector to be populated with the Jacobian entries

Aliasing restrictions: res_0 and res cannot alias (obviously).

NonlinearSolvers.calcJacCol

This function extracts the entries for one column of the Jacobian from a complex step residual evaluation

Inputs: res: vector of perturbed residual values epsilon: magnitude of perturbation

Inputs/Outputs: jac_row = vector to be populated with the Jacobian entries

Aliasing restrictions: none

Explicit Calculation

The most efficient method to compute the Jacobian is to explicitly compute all its entries and assemble them into a matrix. The evalJacobian function must be extended by each physics module for this to work. This function is passed all the same arguments as evalResidual plus an additional _AssembleElementData object which is used to assemble the contribution of each element or interface into the matrix.

NonlinearSolvers._AssembleElementData — Type.

Helper object for assembling element and interface jacobians into the system matrix.

Fields

A: the matrix (can be Array, SparseMatrixCSC, or PetscMat)
idx: temporary array for row indices, length numDofPerNode
idy: temporary array for column indices, length numDofPerNode
vals: temporary array for matrix entries, size numDofPerNode square
idx_i: temporary array for row indices when assembling interface jacobians, length 2 x numDofPerNode
idy_i: like idx_i, but for column indices
vals_i: temporary array for storing matrix entries when assembling interface jacobians, size 2 x numDofPerNode square

NonlinearSolvers._AssembleElementData — Method.

Outer constructor for _AssembleElementData

NonlinearSolvers.NullAssembleElementData — Constant.

An empty _AssembleElementData. Useful for giving a default value to fields.

NonlinearSolvers.assembleElement — Method.

Inputs

helper: an _AssembleElementData
mesh: a mesh
elnum: element number
jac: 4 dimensional array containing the jacobian of the element

jac contains the data for the jacobian of the volume terms for a given element. $jac[i, j, p, q] = \partial R[i, p, elnum] / \partial eqn.q[j, q, elnum]$. Its size is numDofPerNode x numDofPerNode x numNodesPerElement x numNodesPerElement.

NonlinearSolvers.assembleInterface — Function.

Assembles contributions into an element-block matrix. jacLR and jacRL are not used.

Assembles the jacobian of an interface into the matrix. Specialized versions take advantage of the sparsity of the sbpface.

Inputs

helper: _AssembleElementData
sbpface: an SBP face object
mesh: a mesh
iface: an Interface object identify the interface to be assembled
jacLL: see below
jacLR:
jacRL
jacRR

jacAB where A = L or R and B = L or R, is the jacobian of the residual of element A with respect to the solution of element B.

jacAB has the same size/layout as jac in assembleElement.

NonlinearSolvers.assembleSharedFace — Function.

Assembles contribution of the shared face terms into the element-block matrix. jacLR is not used.

Assemble one half of an interface, used by shared face integrals. See assembleInterface.

Inputs

helper: _AssembleElementData
sbpface: an SBP face object
mesh: a mesh
jacLL
jacLR

NonlinearSolvers.assembleBoundary — Function.

Assembles the boundary terms into an element-block matrix.

Assembles the jacobian of a boundary integral into the matrix. Specialized versions take advantage of the sparsity of the sbpface.

Inputs

helper: _AssembleElementData
sbpface: an SBP face object
mesh: a mesh
jac: a jac, same layout as assembleElement