Numerical Integration

  • The only remaining non-discretized parts of the weak form are the integrals.

  • We split the domain integral into a sum of integrals over elements:

  • Through a change of variables, the element integrals are mapped to integrals over the "reference" elements .

  • is the Jacobian of the map from the physical element to the reference element.

  • To approximate the reference element integrals numerically, we use quadrature (typically "Gaussian Quadrature"):

  • is the spatial location of the th quadrature point and is its associated weight.

  • MOOSE handles multiplication by the Jacobian and the weight automatically, thus your Kernel is only responsible for computing the part of the integrand.

  • Under certain common situations, the quadrature approximation is exact! - For example, in 1 dimension, Gaussian Quadrature can exactly integrate polynomials of order with quadrature points.

  • Note that sampling at the quadrature points yields:

  • And our weak form becomes:

  • The second sum is over boundary faces, .

  • MOOSE Kernels must provide each of the terms in square brackets (evaluated at or as necessary).

Newton's Method

  • We now have a nonlinear system of equations, to solve for the coefficients .

  • Newton's method has good convergence properties, we use it to solve this system of nonlinear equations.

  • Newton's method is a "root finding" method: it finds zeros of nonlinear equations.

  • Newton's Method in "Update Form" for finding roots of the scalar equation

  • We don't have just one scalar equation: we have a system of nonlinear equations.

  • This leads to the following form of Newton's Method:

  • Where is the Jacobian matrix evaluated at the current iterate:

  • Note that:

Newton for a Simple Equation

  • Consider the convection-diffusion equation with nonlinear , , and :

  • The component of the residual vector is:

  • Using the previously-defined rules for and , the entry of the Jacobian is then:

  • Note that even for this "simple" equation, the Jacobian entries are nontrivial: they depend on the partial derivatives of , , and , which may be difficult or time-consuming to compute analytically.

  • In a multiphysics setting with many coupled equations and complicated material properties, the Jacobian might be extremely difficult to determine.

Chain Rule

  • On the previous slide, the term was used, where was a nonlinear forcing function.

  • The chain rule allows us to write this term as

  • If a functional form of is known, e.g. , this formula implies that its Jacobian contribution is given by

Jacobian Free Newton Krylov

  • is a linear system solved during each Newton step.

  • For simplicity, we can write this linear system as , where: - - -

  • We employ an iterative Krylov method (e.g. GMRES) to produce a sequence of iterates ,

  • and remain fixed during the iterative process.

  • The "linear residual" at step is defined as

  • MOOSE prints the norm of this vector, , at each iteration, if you set print_linear_residuals = true in the Outputs block.

  • The "nonlinear residual" printed by MOOSE is .

  • By iterate , the Krylov method has constructed the subspace

  • Different Krylov methods produce the iterates in different ways:

  • Conjugate Gradients: orthogonal to .

  • GMRES/MINRES: has minimum norm for in .

  • Biconjugate Gradients: is orthogonal to

  • is never explicitly needed to construct the subspace, only the action of on a vector is required.

  • This action can be approximated by:

  • This form has many advantages: - No need to do analytic derivatives to form - No time needed to compute (just residual computations) - No space needed to store

Wrap Up

  • The Finite Element Method is a way of numerically approximating the solution of PDEs.

  • Just like polynomial fitting, FEM finds coefficients for basis functions.

  • The "solution" is the combination of the coefficients and the basis functions, and the solution can be sampled anywhere in the domain.

  • We compute integrals numerically using quadrature.

  • Newton's Method provides a mechanism for solving a system of nonlinear equations.

  • The Jacobian Free Newton Krylov (JFNK) method allows us to avoid explicitly forming the Jacobian matrix while still computing its "action".