/[escript]/trunk/doc/user/symbolic.tex
ViewVC logotype

Contents of /trunk/doc/user/symbolic.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4888 - (show annotations)
Mon Apr 28 00:52:20 2014 UTC (5 years, 2 months ago) by jfenwick
File MIME type: application/x-tex
File size: 16130 byte(s)
Removing manual pagebreaks
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(yconstraint-1),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}{|c|c|c|}
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 Newton-Raphson 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 Newton-Raphson 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}

  ViewVC Help
Powered by ViewVC 1.1.26