From 4376add9ec2b9fc326e1355645c1208a174285c0 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 31 Jul 2018 15:55:53 +0200
Subject: [PATCH] Laplace Eq. test and docu
---
.../bcs_laplace_eq.py | 52 ++++++++
.../python_laplace_eq_ref.vtu | 3 +
.../square_1e3_laplace_eq.prj | 122 ++++++++++++++++++
.../square_1x1.gml | 3 +
.../square_1x1_quad_1e3.vtu | 3 +
.../laplace-equation/python-laplace-eq.md | 95 ++++++++++++++
.../python_laplace_eq_diff.png | 3 +
.../python_laplace_eq_solution.png | 3 +
8 files changed, 284 insertions(+)
create mode 100644 Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py
create mode 100644 Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu
create mode 100644 Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj
create mode 100644 Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
create mode 100644 Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
create mode 100644 web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md
create mode 100644 web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png
create mode 100644 web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py
new file mode 100644
index 00000000000..e6d5c8dd7d8
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/bcs_laplace_eq.py
@@ -0,0 +1,52 @@
+
+import OpenGeoSys
+from math import pi, sin, cos, sinh, cosh
+
+a = 2.0*pi/3.0
+
+# analytical solution used to set the Dirichlet BCs
+def solution(x, y):
+ return sin(a*x) * sinh(a*y)
+
+# gradient of the analytical solution used to set the Neumann BCs
+def grad_solution(x, y):
+ return a * cos(a*x) * sinh(a*y), \
+ a * sin(a*x) * cosh(a*y)
+
+# Dirichlet BCs
+class BCTop(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert y == 1.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+class BCLeft(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert x == 0.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+class BCBottom(OpenGeoSys.BoundaryCondition):
+ def getDirichletBCValue(self, t, coords, node_id, primary_vars):
+ x, y, z = coords
+ assert y == 0.0 and z == 0.0
+ value = solution(x, y)
+ return (True, value)
+
+# Neumann BC
+class BCRight(OpenGeoSys.BoundaryCondition):
+ def getFlux(self, t, coords, primary_vars):
+ x, y, z = coords
+ assert x == 1.0 and z == 0.0
+ value = grad_solution(x, y)[0]
+ Jac = [ 0.0 ] # value does not depend on primary variable
+ return (True, value, Jac)
+
+
+# instantiate BC objects referenced in OpenGeoSys' prj file
+bc_top = BCTop()
+bc_right = BCRight()
+bc_bottom = BCBottom()
+bc_left = BCLeft()
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu
new file mode 100644
index 00000000000..ca30c517242
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/python_laplace_eq_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5afedc3c451f400d58bdc8cefe7d5989bc1a248d80470131f750b2eb4793e49f
+size 23684
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj
new file mode 100644
index 00000000000..5e196baa339
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj
@@ -0,0 +1,122 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>square_1x1_quad_1e3.vtu</mesh>
+ <geometry>square_1x1.gml</geometry>
+ <python_script>bcs_laplace_eq.py</python_script>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+
+ <jacobian_assembler>
+ <type>CentralDifferences</type>
+ </jacobian_assembler>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_newton</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>square_1e3_neumann</prefix>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>zero</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>zero</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Python</type>
+ <bc_object>bc_left</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Python</type>
+ <bc_object>bc_right</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Python</type>
+ <bc_object>bc_top</bc_object>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>square_1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Python</type>
+ <bc_object>bc_bottom</bc_object>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_newton</name>
+ <type>Newton</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
new file mode 100644
index 00000000000..3e5bb2a80a8
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a3f52cd5a8058d967f8ead873e1afa323af39266e538bcd2c35af08683a81dad
+size 1083
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
new file mode 100644
index 00000000000..b56982101f6
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3653db864dd917e8d5be84d68ee2de13ba3485b1ef23de4e9894fb6749a1aa11
+size 25692
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md b/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md
new file mode 100644
index 00000000000..0e68c1d6c6a
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python-laplace-eq.md
@@ -0,0 +1,95 @@
++++
+project = "https://github.com/ufz/ogs-data/blob/master/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_laplace_eq.prj"
+author = "Christoph Lehmann"
+date = "2018-06-01T14:16:55+02:00"
+title = "Manufactured Solution for Laplace's Equation with Python"
+weight = 1
+
+[menu]
+ [menu.benchmarks]
+ parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## 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$.
+
+
+## Results
+
+The numerical result obtained from OpenGeoSys is:
+
+{{< img src="../python_laplace_eq_solution.png" >}}
+
+The absolute difference between the analytical and numerical solutions is
+smaller than $4 \cdot 10^{-4}$:
+
+{{< img src="../python_laplace_eq_diff.png" >}}
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png
new file mode 100644
index 00000000000..5c4d28f5451
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_diff.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:37b86d2f172fa864d6bfecabf99bd203ddaded8210808956325b4a97b3020f25
+size 34437
diff --git a/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png
new file mode 100644
index 00000000000..ad1bab5c5da
--- /dev/null
+++ b/web/content/docs/benchmarks/python-bc/laplace-equation/python_laplace_eq_solution.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d6aa873d3c8cb0f6b4e7cbd30cf32c07d1ae261e313d503a15d6e505902626fc
+size 39593
--
GitLab