Project file on GitHub
## Poisson Equation

## Problem specification and analytical solution

## Input files

## Running simulation

## Results and evaluation

### Comparison of the numerical and analytical solutions

#### Comparison for higher resolution mesh (\(316 \times 316\) elements)

The Poisson equation is: \[ \begin{equation} - k\; \Delta p = Q \quad \text{in }\Omega \end{equation}\] w.r.t boundary conditions \[ \eqalign{ p(x) = g_D(x) &\quad \text{on }\Gamma_D,\cr k\;{\partial p(x) \over \partial n} = g_N(x) &\quad \text{on }\Gamma_N, }\]

where \(p\) could be the pressure, the subscripts \(D\) and \(N\) denote the Dirichlet- and Neumann-type boundary conditions, \(n\) is the normal vector pointing outside of \(\Omega\), and \(\Gamma = \Gamma_D \cup \Gamma_N\) and \(\Gamma_D \cap \Gamma_N = \emptyset\).

We solve the Poisson equation on a square domain \(\Omega = [0\times 1]^2\) with \(k = 1\) w.r.t. the specific boundary conditions: \[ \eqalign{ p(x,y) = 1 &\quad \text{on } (x=0,y=0) \subset \Gamma_D,\cr p(x,y) = 1 &\quad \text{on } (x=0,y=1) \subset \Gamma_D,\cr p(x,y) = 1 &\quad \text{on } (x=1,y=0) \subset \Gamma_D,\cr p(x,y) = 1 &\quad \text{on } (x=1,y=1) \subset \Gamma_D,\cr }\] and the source term is \[ Q(x,y) = a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) +b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) \]

The analytical solution of (1) is \[ p(x,y) = \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right). \] Since \[ \frac{\partial^2 p}{\partial x} (x,y) = - a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) \] and \[ \frac{\partial^2 p}{\partial y} (x,y) = - b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) \] it yields \[ \Delta p(x,y) = - a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) - b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) \] and \[ Q = - \Delta p(x,y) = a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right) + b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right). \]

The project file is `square_1e3_poisson_sin_x_sin_y.prj`

. It describes the process to be solved and the related process variable together with their initial, boundary conditions and source terms. The definition of the source term \(Q\) is in a Python script. The script for setting the source terms is referenced in the project file as follows:

`<python_script>sin_x_sin_y_source_term.py</python_script>`

In the source term descripition

```
<source_term>
<mesh>square_1x1_quad_1e3_entire_domain</mesh>
<type>Python</type>
<source_term_object>sinx_siny_source_term</source_term_object>
</source_term>
```

the domain is specified by the mesh-tag. The function \(Q\) is defined by the Python object `sinx_sinx_source_term`

that is created in the last line of the Python script:

```
import OpenGeoSys
from math import pi, sin
a = 2.0*pi
b = 2.0*pi
def solution(x, y):
return - sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
# - laplace(solution) = source term
def laplace_solution(x, y):
return a*a * sin(a*x-pi/2.0) * sin(b*y-pi/2.0) + b*b * sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
# source term for the benchmark
class SinXSinYSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
value = laplace_solution(x,y)
Jac = [ 0.0 ]
return (value, Jac)
# instantiate source term object referenced in OpenGeoSys' prj file
sinx_siny_source_term = SinXSinYSourceTerm()
```

To start the simulation (after successful compilation) run:

It will produce some output and write the computed result into the `square_1e3_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu`

, which can be directly visualized and analysed in paraview for example.

The output on the console will be similar to the following on:

```
info: This is OpenGeoSys-6 version 6.1.0-1132-g00a6062a2.
info: OGS started on 2018-10-10 09:22:17+0200.
info: ConstantParameter: K
info: ConstantParameter: p0
info: ConstantParameter: pressure_edge_points
info: Initialize processes.
createPythonSourceTerm
info: Solve processes.
info: [time] Output of timestep 0 took 0.000695229 s.
info: === Time stepping at step #1 and time 1 with step size 1
info: [time] Assembly took 0.0100119 s.
info: [time] Applying Dirichlet BCs took 0.000133991 s.
info: ------------------------------------------------------------------
info: *** Eigen solver computation
info: -> solve with CG (precon DIAGONAL)
info: iteration: 81/10000
info: residual: 5.974447e-17
info: ------------------------------------------------------------------
info: [time] Linear solver took 0.00145817 s.
info: [time] Iteration #1 took 0.0116439 s.
info: [time] Solving process #0 took 0.011662 s in time step #1
info: [time] Time step #1 took 0.0116858 s.
info: [time] Output of timestep 1 took 0.000671864 s.
info: The whole computation of the time stepping took 1 steps, in which
the accepted steps are 1, and the rejected steps are 0.
info: [time] Execution took 0.0370049 s.
info: OGS terminated on 2018-10-10 09:22:17+020
```

The above picture shows the numerical result. The solution conforms in the edges to the prescribed boundary conditions. Since a coarse mesh (\(32 \times 32\) elements) is used for the simulation the difference between the numerical and the analytical solution is relatively large.

The difference between the numerical and the analytical solution is much smaller than in the coarse mesh case.

Previous

This article was written by Tom Fischer. If you are missing something or you find an error please let us know.
Generated with Hugo 0.54.0. Last revision: December 3, 2018 *
Commit: [web] Fixed links. 651493dbb
| Edit this page on
*