1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2003-2018 by The University of Queensland
% http://www.uq.edu.au
%
% Primary Business: Queensland, Australia
% Licensed under the Apache License, version 2.0
% http://www.apache.org/licenses/LICENSE-2.0
%
% Development until 2012 by Earth Systems Science Computational Center (ESSCC)
% Development 2012-2013 by School of Earth Sciences
% Development from 2014 by Centre for Geoscience Computing (GeoComp)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{The \ripley Module}\label{chap:ripley}\index{ripley}
%\declaremodule{extension}{ripley}
%\modulesynopsis{Solving linear, steady partial differential equations using finite elements}
\ripley is an alternative domain library to \finley; it supports structured, uniform meshes with rectangular elements in 2D and hexahedral elements in 3D.
Uniform meshes allow a straightforward division of elements among processes
with \MPI and allow for a number of optimizations when solving PDEs.
\ripley also supports fast assemblers for certain types of PDE (specifically
Lam\'e and Wave PDEs).
These assemblers make use of the regular nature of the domain to optimize the
stiffness matrix assembly process for these specific problems.
Finally, \ripley is currently the only domain family that supports GPU-based
solvers.
As a result, \ripley domains cannot be created by reading from a mesh file
since only one element type is supported and all elements need to be equally
sized.
For the same reasons, \ripley does not allow assigning coordinates via
\function{setX}.
While \ripley cannot be used with mesh files, it can be used to read in \GOCAD \index{GOCAD} \index{ripley!GOCAD}
data. A script with an example of a voxet reader is included in the examples as
\texttt{voxet_reader.py}.
Other than use of meshfiles, \ripley and \finley are generally interchangeable in a script
with both modules having the \class{Rectangle} or \class{Brick} functions
available. Consider the following example which creates a 2D \ripley domain:
\begin{python}
from esys.ripley import Rectangle, Brick
dom = Rectangle(9, 9)
\end{python}
\index{ripley!multi-resolution domains} \index{multi-resolution domains}
Multi-resolution domains are supported in \ripley via \class{MultiBrick}
and \class{MultiRectangle}. Each level of one of
these domains has twice the elements in each axis of the next lower resolution.
The \class{MultiBrick} is not currently supported when running \escript with
multiple processes using \MPI. Interpolation between these multi-resolution
domains is possible providing they have matching dimensions and subdivisions,
along with a compatible number of elements.
To simplify these conditions the
use of \class{MultiResolutionDomain} is highly recommended. The following
example creates two 2D domains of different resolutions and interpolates
between them:
\begin{python}
from esys.ripley import MultiResolutionDomain
mrd = MultiResolutionDomain(2, n0=10, n1=10)
ten_by_ten = mrd.getLevel(0)
data10 = Vector(..., Function(ten_by_ten))
...
forty_by_forty = mrd.getLevel(2)
data40 = interpolate(data10, Function(forty_by_forty))
\end{python}
\section{Formulation}
For a single PDE that has a solution with a single component the linear PDE is
defined in the following form:
\begin{equation}\label{eq:ripleysingle}
\begin{array}{cl} &
\displaystyle{
\int_{\Omega}
A_{jl} \cdot v_{,j}u_{,l}+ B_{j} \cdot v_{,j} u+ C_{l} \cdot v u_{,l}+D \cdot vu \; d\Omega }
+ \displaystyle{\int_{\Gamma} d \cdot vu \; d{\Gamma} }\\
= & \displaystyle{\int_{\Omega} X_{j} \cdot v_{,j}+ Y \cdot v \; d\Omega }
+ \displaystyle{\int_{\Gamma} y \cdot v \; d{\Gamma}}
\end{array}
\end{equation}
\section{Meshes}
\label{sec:ripleymeshes}
An example 2D mesh from \ripley is shown in Figure~\ref{fig:ripleyrect}.
Mesh files cannot be used to generate \ripley domains, i.e. \ripley does not
have \function{ReadGmsh} or \function{ReadMesh} functions.
Instead, \ripley domains are always created using a call to \function{Brick}
or \function{Rectangle}, see Section~\ref{sec:ripleyfuncs}.
\begin{figure}
\centerline{\includegraphics{RipleyMesh}}
\caption{10x10 \ripley Rectangle created with \var{l0=(5.5, 15.5)} and
\var{l1=(9.0, 14.0)}}
\label{fig:ripleyrect}
\end{figure}
\section{Functions}\label{sec:ripleyfuncs}
\begin{funcdesc}{Brick}{n0,n1,n2,l0=1.,l1=1.,l2=1.,d0=-1,d1=-1,d2=-1,
diracPoints=list(), diracTags=list()}
generates a \Domain object representing a three-dimensional brick between
$(0,0,0)$ and $(l0,l1,l2)$ with orthogonal faces. All elements will be regular.
The brick is filled with
\var{n0} elements along the $x_0$-axis,
\var{n1} elements along the $x_1$-axis and
\var{n2} elements along the $x_2$-axis.
If built with \MPI support, the domain will be subdivided
\var{d0} times along the $x_0$-axis,
\var{d1} times along the $x_1$-axis, and
\var{d2} times along the $x_2$-axis.
\var{d0}, \var{d1}, and \var{d2} must be factors of the number of
\MPI processes requested.
If axial subdivisions are not specified, automatic domain subdivision will take
place. This may not be the most efficient construction and will likely result in
extra elements being added to ensure proper distribution of work. Any extra
elements added in this way will change the length of the domain proportionately.
\var{diracPoints} is a list of coordinate-tuples of points within the mesh,
each point tagged with the respective string within \var{diracTags}.
\end{funcdesc}
\begin{funcdesc}{Rectangle}{n0,n1,l0=1.,l1=1.,d0=-1,d1=-1,
diracPoints=list(), diracTags=list()}
generates a \Domain object representing a two-dimensional rectangle between
$(0,0)$ and $(l0,l1)$ with orthogonal faces. All elements will be regular.
The rectangle is filled with
\var{n0} elements along the $x_0$-axis and
\var{n1} elements along the $x_1$-axis.
If built with \MPI support, the domain will be subdivided
\var{d0} times along the $x_0$-axis and
\var{d1} times along the $x_1$-axis.
\var{d0} and \var{d1} must be factors of the number of \MPI processes requested.
If axial subdivisions are not specified, automatic domain subdivision will take
place. This may not be the most efficient construction and will likely result in
extra elements being added to ensure proper distribution of work. Any extra
elements added in this way will change the length of the domain proportionately.
\var{diracPoints} is a list of coordinate-tuples of points within the mesh,
each point tagged with the respective string within \var{diracTags}.
\end{funcdesc}
\noindent The arguments \var{l0}, \var{l1} and \var{l2} for \function{Brick}
and \function{Rectangle} may also be given as tuples \var{(x0,x1)} in which
case the coordinates will range between \var{x0} and \var{x1}. For example:
\begin{python}
from esys.ripley import Rectangle
dom = Rectangle(10, 10, l0=(5.5, 15.5), l1=(9.0, 14.0))
\end{python}
\noindent This will create a rectangle with $10$ by $10$ elements where the
bottom-left node is located at $(5.5, 9.0)$ and the top-right node has
coordinates $(15.5, 14.0)$, see Figure~\ref{fig:ripleyrect}.
The \class{MultiResolutionDomain} class is available as a wrapper, taking the
dimension of the domain followed by the same arguments as \class{Brick} (if
a two-dimensional domain is requested, any extra arguments over those used by
\class{Rectangle} are ignored). All of these standard arguments to
\class{MultiResolutionDomain} must be supplied as keyword arguments
(e.g. \var{d0}=...). The \class{MultiResolutionDomain} can then generate
compatible domains for interpolation.
\section{Linear Solvers in \SolverOptions}
Currently direct solvers and GPU-based solvers are not supported under \MPI
when running with more than one rank.
By default, \ripley uses the iterative solvers \PCG for symmetric and \BiCGStab
for non-symmetric problems.
A GPU will not be used unless explicitly requested via the
\function{setSolverTarget} method of the solver options.
These solvers are only available if \ripley was built with \CUDA support.
If the direct solver is selected, which can be useful when solving very
ill-posed equations, \ripley uses the \MKL\footnote{If the stiffness matrix is
non-regular \MKL may return without a proper error code. If you observe
suspicious solutions when using \MKL, this may be caused by a non-invertible
operator.} solver package. If \MKL is not available \UMFPACK is used.
If \UMFPACK is not available a suitable iterative solver from \PASO is used, but
if a direct solver was requested via \SolverOptions an exception will be
raised.
|
