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 Neumanntype 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 description
<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 meshtag. 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*xpi/2.0) * sin(b*ypi/2.0)
#  laplace(solution) = source term
def laplace_solution(x, y):
return a*a * sin(a*xpi/2.0) * sin(b*ypi/2.0) + b*b * sin(a*xpi/2.0) * sin(b*ypi/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:
ogs square_1e3_poisson_sin_x_sin_y.prj
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 OpenGeoSys6 version 6.1.01132g00a6062a2.
info: OGS started on 20181010 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.974447e17
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 20181010 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.
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.101.0
in CI job 319943

Last revision: February 16, 2023
Commit: [PL] Extract another loop and collection of IVs de82977c1
 Edit this page on