1 
%!TEX root = user.tex 
2 
\chapter{Symbolic} 
3 
\label{CHAP: The escript symbolic toolbox} 
4 
\section{Introduction} 
5 
Escript builds on the existing Sympy\cite{Sympy} symbolic maths library to provide a \SYMBOL class with support for escript Data objects. \SYMBOL objects act as place holder's for a single mathematical symbol, 
6 
such as x, or for arbitrarily complex mathematical expressions such as 
7 
c * x**4 + alpha * exp(x)  2 * sin(beta * x), where alpha, beta, c, and x 
8 
are also Symbols (the symbolic ``atoms" of the expression). 
9 
With the help of the \EVALUATOR class, these symbols and expressions can 
10 
be evaluated by substituting numeric values and/or escript \Data objects 
11 
for the atoms. Escripts \SYMBOL class has a 
12 
shape (and thus a rank) as well as a dimension. 
13 
Symbols are useful to perform mathematical simplifications, 
14 
compute derivatives, take gradients and in the case of escript describe PDE's. As an example of how the symbolic toolbox can be used, consider the following code extract. 
15 
\begin{python} 
16 
import esys.escript as es 
17 
u = es.Symbol('u') 
18 
p = 2*u**2+3*u+1 
19 
p2 = es.sin(u) 
20 
p3 = p.diff(u) #p3 = derivative of p with respect to u 
21 
evalu = es.Evaluator() 
22 
evalu.addExpression(p) 
23 
evalu.addExpression(p2) 
24 
evalu.addExpression(p3) 
25 
evalu.subs(u=2*es.symconstants.pi) 
26 
evaluated=evalu.evaluate() 
27 
print evaluated 
28 
\end{python} 
29 
Running this code outputs (1 + 6*pi + 8*pi**2, 0, 3 + 8*pi). To get the numeric value of the expression we replace evalu.evaluate() with evalu.evaluate(evalf=True) this results in (98.8063911302536, 0, 28.1327412287183). The use of these Symbols becomes more interesting in the context of escript when they are integrated with escript \Data objects. 
30 
\section{NonlinearPDE} 
31 
The NonlinearPDE class in escript makes use of the escript \SYMBOL class. The NonlinearPDE class allows for the solution of PDEs of the form: 
32 
\begin{equation} 
33 
div(X) + Y = 0 
34 
\label{symbolic eq1} 
35 
\end{equation} 
36 
where $X$ and $Y$ are both functions of $grad(u)$ and $u$. Where $u$ is the unknown function implemented as a \SYMBOL. div(x) denotes divergence of x. 
37 
The NonlinearPDE class uses the Symbolic class to solve the nonlinear PDE given in \eqn{symbolic eq1} 
38 
The class works by using Newton's method to find the zeroes of the left hand side of \eqn{symbolic eq1} and as a consequence finding the 
39 
$X$ and $Y$ which satisfy \eqn{symbolic eq1}. 
40 
Consecutive updates are calculated until the equation is 
41 
satisfied to the desired level of accuracy. The solution to each update step involves solving a linear PDE. The NonlinearPDE class uses $X$ and $Y$ to produce the coefficient of the linear PDE for the update step. The linear PDE class given in \Sec{SEC LinearPDE} is used to solve the linear PDEs from the update step. The coefficients of the linear PDE to be solved are calculated as follows: 
42 
\begin{equation*} 
43 
A = \frac{\partial \text{X}}{\partial grad(u)}, B = \frac{\partial \text{X}}{\partial u}, C = \frac{\partial \text{Y}}{\partial grad(u)}, D = \frac{\partial \text{Y}}{\partial u} 
44 
\end{equation*} 
45 
\section{2D Plane Strain Problem} 
46 
The NonLinearPDE class can be used to solve a 2D plane strain problem. In continuous media, stress is given by Lame's equation, \eqn{symbolic eq2}. 
47 
\begin{equation} 
48 
div(\sigma)=f 
49 
\label{symbolic eq2} 
50 
\end{equation} 
51 
Hook's Law provides a relation between $\sigma$ and $\epsilon$ in the following form 
52 
\begin{equation} 
53 
\left[ \begin{array}{c} 
54 
\sigma_{00} \\ 
55 
\sigma_{11} \\ 
56 
\sigma_{01} \\ 
57 
\end{array} \right] = 
58 
\left[ \begin{array}{ccc} 
59 
c_{00} & c_{01} & c_{05}\\ 
60 
c_{01} & c_{11} & c_{15}\\ 
61 
c_{05} & c_{15} & c_{55}\\ 
62 
\end{array}\right] 
63 
\left[ \begin{array}{c} 
64 
\epsilon_{00} \\ 
65 
\epsilon_{11} \\ 
66 
2\epsilon_{10} \\ 
67 
\end{array} \right] 
68 
\label{symbolic eq3} 
69 
\end{equation} 
70 
Where $\epsilon = symmetric(grad(u)) \text{ or } \epsilon_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + {\frac{\partial u_j}{\partial x_i}}\right)$, u is the unknown function and $c_{ij}$ is the stiffness matrix. To fit this to the nonlinear PDE class' standard form, X, is set to $c \times Symmetric(\text{grad(}u)))$. 
71 
The following python extract show how an example 2D plane strain problem can be set up. 
72 

73 

74 
\begin{python} 
75 
from esys.escript import * 
76 
from esys.finley import Rectangle 
77 
#set up domain and symbols 
78 
mydomain = Rectangle(l0=1.,l1=1.,n0=10, n1=10) 
79 
u = Symbol('u',(2,), dim=2) 
80 
q = Symbol('q', (2,2)) 
81 
sigma = Symbol('sigma',(2,2)) 
82 
theta = Symbol('theta') 
83 
# q is a rotation matrix represented by a Symbol. Values can be substituted for 
84 
# theta. 
85 
q[0,0]=cos(theta) 
86 
q[0,1]=sin(theta) 
87 
q[1,0]=sin(theta) 
88 
q[1,1]=cos(theta) 
89 
# Theta gets substituted by pi/4 and masked to lie between .3 and .7 in the 
90 
# vertical direction. Using this masking means that when q is used it will apply 
91 
# only to the specified area of the domain. 
92 
x = Function(mydomain).getX() 
93 
q=q.subs(theta,(symconstants.pi/4)*whereNonNegative(x[1].30)*whereNegative(x[1].70)) 
94 
# epsilon is defined in terms of u and has the rotation applied. 
95 
epsilon0 = symmetric(grad(u)) 
96 
epsilon = matrixmult(matrixmult(q,epsilon0),q.transpose(1)) 
97 
# For the purposes of demonstration, an arbitrary c with isotropic constraints 
98 
# is chosen here. In order to act as an isotropic material c is chosen such that 
99 
# c00 = c11 = c01+c1+2*c55 
100 
c00 = 10 
101 
c01 = 8; c11 = 10 
102 
c05 = 0; c15 = 0; c55 = 1 
103 
# sigma is defined in terms of epsilon 
104 
sigma[0,0] = c00*epsilon[0,0]+c01*epsilon[1,1]+c05*2*epsilon[1,0] 
105 
sigma[1,1] = c01*epsilon[0,0]+c11*epsilon[1,1]+c15*2*epsilon[1,0] 
106 
sigma[0,1] = c05*epsilon[0,0]+c15*epsilon[1,1]+c55*2*epsilon[1,0] 
107 
sigma[1,0] = sigma[0,1] 
108 
sigma0=matrixmult(matrixmult(q.transpose(1),sigma),q) 
109 
# set up boundary conditions 
110 
x=mydomain.getX() 
111 
gammaD=whereZero(x[1])*[1,1] 
112 
yconstraint = FunctionOnBoundary(mydomain).getX()[1] 
113 
# The nonlinear PDE is set up, the values are substituted in and the solution is 
114 
# calculated, y represents an external shearing force acting on the domain. 
115 
# In this case a force of magnitude 50 is acting in the x[0] direction. 
116 
p = NonlinearPDE(mydomain, u, debug=NonlinearPDE.DEBUG0) 
117 
p.setValue(X=sigma0,q=gammaD,y=[50,0]*whereZero(yconstraint1),r=[1,1]) 
118 
v = p.getSolution(u=[0,0]) 
119 
\end{python} 
120 
%\pagebreak 
121 
The way in which the rotation matrix q is set up demonstrates the seamless integration of escript symbols and \Data objects. A \SYMBOL is used to set up the matrix, the values for theta are then later substituted in. The example also demonstrates how the symbolic toolbox can be used as an aid to easily move from a mathematical equation to an escript data object which can be used to do numerical calculations. 
122 
Running the script calculates the unknown function u and assigns it to v. We can use v to calculate the stress and strain. 
123 
\begin{table}[!h] 
124 
\centering 
125 
\begin{tabular}{ccc} 
126 
\hline 
127 
& Anisotropic & Isotropic\\%\multicolumn{2}{c}{Isotopic} \\ 
128 
\hline 
129 
No rotation & \includegraphics[scale=0.2]{0RotAniso} & \includegraphics[scale=0.2]{0RotIso}\\ 
130 
\hline 
131 
60\textdegree\ rotation & \includegraphics[scale=0.2]{MidRotAniso} & \includegraphics[scale=0.2]{MidRotIso}\\ 
132 
\hline 
133 
difference & \includegraphics[scale=0.2]{diffaniso} & \raisebox{2cm}{no difference}\\ 
134 
\hline 
135 
\end{tabular} 
136 
\caption{Displacement vectors calculated using NonlinearPDE} 
137 
\label{isovsaniso} 
138 
\end{table} 
139 
Table \ref{isovsaniso}, shows the result of running the above script under varying values of c and theta. Both isotropic and anisotropic cases are considered. For the anisotropic case, $c_{ij}$ is chosen such that c00 = c11 = c01+c1+2*c55 does not hold. Two values of theta are also considered; one with a masked 60\textdegree rotation in the middle and one with no rotation. The last row of the table shows the difference between rotation in the middle and no rotation. In the isotropic case it can be seen that there is no difference in the output when the rotation is applied. There is however, an obvious difference when the anisotropic case is considered. 
140 

141 
\newpage 
142 
\section{Classes} 
143 
A number of classes are associated with escript symbols. A detailed listing of the definitions and usage is provided below. 
144 
\subsection{Symbol class} 
145 
\begin{classdesc}{Symbol}{symbol \optional{, shape} \optional{, Dim}} 
146 
Defines a \SYMBOL object. The first argument \SYMBOL is a string given to represent the \SYMBOL. The string typically matches the name of the object, for instance \texttt{u=Symbol('u')}. Next the shape defines whether the \SYMBOL is a vector or a matrix etc up to a 4 tensor and the length or size of it. Dim is used to define the dimensionality of the object contained in the \SYMBOL. 
147 
For a \SYMBOL definition u = Symbol('u',(10,), dim=2), the value of u will be a vector with 10 components. 
148 
The domain on which u will be used is 2 dimensional(this is relevant with operations such as grad where the number of spacial direction is important). 
149 
\end{classdesc} 
150 
\subsubsection{Symbol class methods} 
151 
\begin{methoddesc}[Symbol]{atoms}{\optional{types}} 
152 
Returns the atoms that form the current \SYMBOL. 
153 
By default, only objects that are truly atomic and cannot be divided 
154 
into smaller pieces are returned: symbols, numbers, and number 
155 
symbols like e and pi. It is possible to request atoms of any type, 
156 
however. 
157 
\end{methoddesc} 
158 

159 
\begin{methoddesc}[Symbol]{coeff}{x \optional{, expand=true}} 
160 
Returns the coefficient of the term x or 0 if there is no x. 
161 
If x is a scalar \SYMBOL then x is searched in all components of 
162 
this \SYMBOL. Otherwise the shapes must match and the coefficients are 
163 
checked component by component. 
164 
Example: 
165 
\begin{python} 
166 
x=Symbol('x', (2,2)) 
167 
y=3*x 
168 
print y.coeff(x) 
169 
print y.coeff(x[1,1]) 
170 
\end{python} 
171 
will print: 
172 
\begin{python} 
173 
[[3 3] 
174 
[3 3]] 
175 
[[0 0] 
176 
[0 3]] 
177 
\end{python} 
178 
\end{methoddesc} 
179 
\begin{methoddesc}[Symbol]{diff}{symbols} 
180 
Takes the derivative of the \SYMBOL object which the method belongs to with respect to the symbols specified in the argument symbols. 
181 
\end{methoddesc} 
182 
\begin{methoddesc}[Symbol]{evalf}{} 
183 
Applies the sympy.evalf operation on all elements of the \SYMBOL which are of type or inherit from sympy.Basic. 
184 
\end{methoddesc} 
185 
\begin{methoddesc}[Symbol]{expand}{} 
186 
Applies the sympy.expand operation on all elements in this \SYMBOL. 
187 
\end{methoddesc} 
188 
%\begin{methoddesc}[Symbol]{getDataSubstitutions}{} 
189 
%end{methoddesc} 
190 
\begin{methoddesc}[Symbol]{getDim}{} 
191 
Returns the symbol's spatial dimensionality, or 1 if undefined. 
192 
\end{methoddesc} 
193 
\begin{methoddesc}[Symbol]{getRank}{} 
194 
Returns the symbol's rank which is equal to the length of the shape. 
195 
\end{methoddesc} 
196 
\begin{methoddesc}[Symbol]{getShape}{} 
197 
Returns the shape of this \SYMBOL. 
198 
\end{methoddesc} 
199 
\begin{methoddesc}[Symbol]{grad}{\optional{where=none}} 
200 
Returns the gradient of the \SYMBOL to which the function is applied. The \SYMBOL must have a dimensionality defined in order for grad to work. As with the normal escript grad a \FunctionSpace can be specified using the where argument. The \FunctionSpace should be wrapped in a \SYMBOL. To do this, set up a \SYMBOL and then use the `subs' function to substitute in the \FunctionSpace. 
201 
\end{methoddesc} 
202 
\begin{methoddesc}[Symbol]{inverse}{} 
203 
Find the inverse of the \SYMBOL to which the function is applied. Inverse is only valid for square rank 2 symbols. 
204 
\end{methoddesc} 
205 
%\begin{methoddesc}[Symbol]{item}{} 
206 
%\end{methoddesc} 
207 
%\begin{methoddesc}[Symbol]{lambdarepr}{} 
208 
%test 
209 
%\end{methoddesc} 
210 
\begin{methoddesc}[Symbol]{simplify}{} 
211 
Applies the sympy.simplify operation on all elements in this \SYMBOL. 
212 
\end{methoddesc} 
213 
\begin{methoddesc}[Symbol]{subs}{old,new} 
214 
Substitutes or replaces a \SYMBOL specified in old with whatever is in new for the \SYMBOL which it is called on. Consider: 
215 
\begin{python} 
216 
import esys.escript as es 
217 
u=es.Symbol("u") 
218 
expr=2*u 
219 
expr.subs(u,2) 
220 
\end{python} 
221 
This prints 4. 
222 
\end{methoddesc} 
223 
\begin{methoddesc}[Symbol]{trace}{axis_offset} 
224 
Returns the trace of the \SYMBOL object. 
225 
\end{methoddesc} 
226 

227 

228 
\subsection{Evaluator class} 
229 
The \EVALUATOR class is intended to have a group of expressions added to it, substitutions can be made across all expressions and the expressions can then all be evaluated. 
230 
\subsubsection{Evaluator class methods} 
231 
\begin{classdesc}{Evaluator}{\optional{expressions}} 
232 
An \EVALUATOR object is initiated via Evaluator() with an optional argument of expressions to store. 
233 
\end{classdesc} 
234 
\begin{methoddesc}[Evaluator]{addExpression}{expression} 
235 
Adds an expression to this \EVALUATOR. 
236 
\end{methoddesc} 
237 
\begin{methoddesc}[Evaluator]{evaluate}{\optional{evalf=False}\optional{, args}} 
238 
Evaluates all expressions in this \EVALUATOR and returns the result as a tuple. Evalf can be set to true to call evalf on any sympy symbols which may be part of the expression. 
239 
Args can be provided to make any substitutions before the expression is evaluated. 
240 
\end{methoddesc} 
241 
\begin{methoddesc}[Evaluator]{subs}{old,new} 
242 
Substitutes or replaces a \SYMBOL specified in old with whatever is in new for all expressions in the \EVALUATOR. 
243 
\end{methoddesc} 
244 

245 
\subsection{NonlinearPDE class} 
246 
\begin{classdesc}{NonlinearPDE}{domain, u} 
247 
Defines a general nonlinear, steady, second order PDE for an unknown function \var{u} on a given domain defined through a \Domain object \var{domain}. \var{u} is a \SYMBOL object. 
248 
The general form is div(X) + Y = 0 
249 
\end{classdesc} 
250 
\iffalse 
251 
\begin{methoddesc}[NonlinearPDE]{concatenateRow}{} 
252 
test 
253 
\end{methoddesc} 
254 
\begin{methoddesc}[NonlinearPDE]{createCoefficient}{} 
255 
test 
256 
\end{methoddesc} 
257 
\begin{methoddesc}[NonlinearPDE]{getUnknownSymbol}{} 
258 
test 
259 
\end{methoddesc} 
260 
\begin{methoddesc}[NonlinearPDE]{getLinearSolverOptions}{} 
261 
test 
262 
\end{methoddesc} 
263 
\begin{methoddesc}[NonlinearPDE]{getLinearPDE}{} 
264 
test 
265 
\end{methoddesc} 
266 
\begin{methoddesc}[NonlinearPDE]{getNumSolutions}{} 
267 
test 
268 
\end{methoddesc} 
269 
\begin{methoddesc}[NonlinearPDE]{getShapeOfCoefficient}{} 
270 
test 
271 
\end{methoddesc} 
272 
\begin{methoddesc}[NonlinearPDE]{getCoefficient}{} 
273 
test 
274 
\end{methoddesc} 
275 
\begin{methoddesc}[NonlinearPDE]{getSensitivity}{} 
276 
test 
277 
\end{methoddesc} 
278 
\fi 
279 
\begin{methoddesc}[NonlinearPDE]{getSolution}{subs} 
280 
Returns the solution of the PDE. Subs contatins the Substitutions for all symbols used in the coefficients including the initial value for the unknown u. 
281 
\end{methoddesc} 
282 
\begin{methoddesc}[NonlinearPDE]{setOptions}{opts} 
283 
\begin{verbatim} 
284 
Allows setting options for the nonlinear PDE. 
285 
The supported options are: 
286 
tolerance 
287 
error tolerance for the Newton method 
288 
iteration_steps_max 
289 
maximum number of Newton iterations 
290 
omega_min 
291 
minimum relaxation factor 
292 
atol 
293 
solution norms less than atol are assumed to be atol. 
294 
This can be useful if one of your solutions is expected to 
295 
be zero. 
296 
quadratic_convergence_limit 
297 
if the norm of the NewtonRaphson correction is reduced by 
298 
less than quadratic_convergence_limit between two iteration 
299 
steps, quadratic convergence is assumed. 
300 
simplified_newton_limit 
301 
if the norm of the defect is reduced by less than 
302 
simplified_newton_limit between two iteration steps and 
303 
quadratic convergence is detected, the iteration switches to the 
304 
simplified NewtonRaphson scheme. 
305 
\end{verbatim} 
306 

307 
\end{methoddesc} 
308 
\begin{methoddesc}[NonlinearPDE]{setValue}{ 
309 
\optional{X}\optional{, Y} 
310 
\optional{, y}\optional{, y_contact} 
311 
\optional{, y_dirac} 
312 
\optional{, q}\optional{, r}} 
313 
Assigns new values to coefficients. By default all values are assumed to be 
314 
zero\footnote{In fact, it is assumed they are not present by assigning the 
315 
value \code{escript.Data()}. This can be used by the solver library to reduce 
316 
computational costs.}. 
317 
If the new coefficient value is not a \Data object, it is converted into a 
318 
\Data object in the appropriate \FunctionSpace. 
319 
\end{methoddesc} 
320 

321 
\subsection{Symconsts class} 
322 
Symconsts provides symbolic constants for use in symbolic expressions. These constants are preferred to floating point implementation as they can cancel perfectly when mathematical expressions are evaluated, avoiding numerical imprecision. 
323 
\begin{verbatim} 
324 
usage: 
325 
symconsts.pi this provides a Symbol object 
326 

327 
Available constants currently include: 
328 
pi and e 
329 
\end{verbatim} 