Manufactured Solution for Laplace's Equation with Python

Project file on GitHub

Motivation of this test case

The aim of this test is:

  • to provide a simple introductory example of how Python BCs can be used
  • to show that Python BCs can be used to easily prescribe BCs based on analytical formulas
  • to check whether both essential and natural BCs are implemented correctly in OpenGeoSys’ Python BC.

Problem description

We solve Laplace’s Equation in 2D on a \(1 \times 1\) square domain. The boundary conditions will be specified below.

Weak form

Laplace’s equation is \[ \begin{equation} - \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) = 0 \end{equation} \] The weak form is derived as usual by multiplying with a test function \(v\) and integrating over the domain \(\Omega\): \[ \begin{equation} - \int_{\Omega} v \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = 0 \,, \end{equation} \] which can be transformed further to \[ \begin{equation} \int_{\Omega} a \mathop{\mathrm{grad}} v \cdot \mathop{\mathrm{grad}} u \, \mathrm{d}\Omega = \int_{\Omega} \mathop{\mathrm{div}} (v a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = \int_{\Gamma_{\mathrm{N}}} v a \mathop{\mathrm{grad}} u \cdot n \, \mathrm{d}\Gamma \,, \end{equation} \] where in the second equality Gauss’s theorem has been applied. As usual, the domain boundary \(\partial\Omega = \Gamma_{\mathrm{D}} \cup \Gamma_{\mathrm{N}}\) is subdivided into the dirichlet and the Neumann boundary and \(v\) vanishes on \(\Gamma_{\mathrm{D}}\). The r.h.s. of the above equation is the total flux associated with \(u\) flowing into the domain \(\Omega\) through \(\Gamma_{\mathrm{N}}\): \(-a \mathop{\mathrm{grad}} u\) is the flux density and \(-n\) is the inwards directed surface normal.

The weak form just derived is implemented (after FEM discretization) in the groundwater flow process in OpenGeoSys. Note that for the application of Neumann boundary conditions, it is necessary to know whether the flux has to be applied with a positive or a negative sign!

Analytial solution

The coefficient \(a\) of Laplace’s equation is taken to be unity. By differentiation it can be easily checked that \[ \begin{equation} u(x, y) = \sin(bx) \sinh(by) \end{equation} \] solves Laplace’s equation inside \(\Omega\) for any \(b\). In this example we set \(b = \tfrac 23 \pi\).

As boundary conditions we apply Dirichlet BCs at the top, left and bottom of the domain with values from \(u(x,y)|_{\Gamma_{\mathrm{D}}}\). On the right boundary of the domain a Neumann BC is applied. There \(n = (1, 0)\), which implies that \(a \mathop{\mathrm{grad}} u \cdot n = a \, \partial u / \partial x\).


The numerical result obtained from OpenGeoSys is:

The absolute difference between the analytical and numerical solutions is smaller than \(4 \cdot 10^{-4}\):

This article was written by Christoph Lehmann. If you are missing something or you find an error please let us know. Generated with Hugo 0.47.1. Last revision: September 12, 2018
Commit: [web] Fix LaTeX formulae with escaped _. 0f264ec51  | Edit this page on

Creative Commons License | Privacy policy / Legal statements / Impressum