Overview

  • A "Kernel" is a piece of physics.
    • It can represent one or more operators or terms in a PDE.
  • A Kernel MUST override computeQpResidual()
  • A Kernel can optionally override:
    • computeQpJacobian()
    • computeQpOffDiagJacobian()

Kernel Member Variables (some)

  • _u, _grad_u
    • Value and gradient of the variable this Kernel is operating on.
  • _phi, _grad_phi
    • Value ($$$\phi$$$) and gradient ($$$\nabla \phi$$$) of the trial functions at the q-points.
  • _test, _grad_test
    • Value ($$$\psi$$$) and gradient ($$$\nabla \psi$$$) of the test functions at the q-points.
  • _q_point
    • XYZ coordinates of the current q-point.
  • _i, _j
    • Current shape functions for test and trial functions, respectively.
  • _qp
    • Current quadrature point index.
  • _current_elem
    • A pointer to the current element that is being operated on.

Diffusion

The strong-form of the diffusion equations is defined as: $$$$\begin{aligned} -\nabla\cdot\nabla u - f &= 0\\ u|_{\partial\Omega_1} &= g_1\\ \nabla u\cdot \hat{n} |_{\partial\Omega_2} &= g_2\end{aligned}$$$$

Generate the weak-form by:

  • Multiplying by test function and integrate over the domain: $$$$(-\nabla\cdot\nabla u, \psi_i) - (f, \psi_i) = 0$$$$
  • Integrating by parts results in the weak form: $$$$(\nabla u, \nabla \psi_i) - (f, \psi_i) - \langle g_2, \psi_i\rangle = 0$$$$

The Jacobian is defined as: $$$$(\nabla \phi_j, \nabla \psi_i)$$$$

Diffusion.h

Diffusion.C

Kernel Registration

  • Before a Kernel is available for use, it must be registered. This is done in e.g. src/base/FrogApp.C for the Frog application.
#include "Diffusion.h"

int MyApp::registerObjects(Factory & factory)
{
  ...
  registerKernel(Diffusion);
  ...
}

Example 2: Adding a Custom Kernel

  • In Example 2, we will register and add a new Convection Kernel
  • Strong Form, excluding BCs: $$$$-\nabla\cdot\nabla u + \vec{v} \cdot\nabla u = 0$$$$
  • Multiply by test function $$$\psi_i$$$ and integrate over the domain: $$$$\left(-\nabla\cdot\nabla u, \psi_i\right) + \left( \vec{v}\cdot\nabla u, \psi_i\right) = 0$$$$
  • Integrate the first term by parts to get our final weak form: $$$$\left(\nabla u, \nabla\psi_i\right) - \langle\nabla u\cdot\hat{n}, \psi_i\rangle + \left(\vec{v} \cdot \nabla u, \psi_i\right) = 0$$$$

Look at Example 2

Time Derivative

The residual contribution for the time derivative term is $$$$\left(\frac{\partial u_h}{\partial t}, \psi_i\right)$$$$ where $$$u_h$$$ is the finite element solution, and $$$$\frac{\partial u_h}{\partial t} \equiv \frac{\partial}{\partial t} \left(\sum_k u_k \phi_k\right) = \sum_k \frac{\partial u_k}{\partial t} \phi_k$$$$ because you can interchange the order of differentiation and summation. Call this equation (1).

In the equation above, $$$\frac{\partial u_k}{\partial t}$$$ is the time derivative of the $$$k$$$th finite element coefficient of $$$u_h$$$. While the exact form of this derivative depends on the time stepping scheme, without much loss of generality, we can assume the following form for the time derivative: $$$$\frac{\partial u_k}{\partial t} = a u_k + b$$$$ for some constants $$$a$$$, $$$b$$$ which depend on $$$\Delta t$$$ and the timestepping method.

The derivative of (1) with respect to $$$u_j$$$ is then:

$$$$\frac{\partial}{\partial u_j} \left(\sum_k \frac{\partial u_k}{\partial t} \phi_k\right) = \frac{\partial }{\partial u_j} \left(\sum_k (a u_k + b) \phi_k\right) = a \phi_j$$$$

So that the jacobian term for (1) is $$$$\left(a \phi_j, \psi_i\right)$$$$ where $$$a$$$ is what we call du_dot_du in MOOSE.

Therefore the computeQpResidual() function in the TimeDerivative Kernel in MOOSE looks like:

return _test[_i][_qp]*_u_dot[_qp];

And the corresponding computeQpJacobian() is:

return _test[_i][_qp]*_phi[_j][_qp]*_du_dot_du[_qp];