ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 4886 - (show annotations)
Sun Apr 27 13:12:23 2014 UTC (5 years, 4 months ago) by jduplessis
File MIME type: application/x-tex
File size: 15989 byte(s)
rewording and checking

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}b$ 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.
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} & 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 ansisotropic 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.
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), u will be a 10 component vector of 2 dimensional objects.
148 \end{classdesc}
149 \subsubsection{Symbol class methods}
150 \begin{methoddesc}[Symbol]{atoms}{\optional{types}}
151 Returns the atoms that form the current \SYMBOL.
152 By default, only objects that are truly atomic and cannot be divided
153 into smaller pieces are returned: symbols, numbers, and number
154 symbols like e and pi. It is possible to request atoms of any type,
155 however.
156 \end{methoddesc}
158 \begin{methoddesc}[Symbol]{coeff}{x \optional{, expand=true}}
159 Returns the coefficient of the term x or 0 if there is no x.
160 If x is a scalar \SYMBOL then x is searched in all components of
161 this \SYMBOL. Otherwise the shapes must match and the coefficients are
162 checked component by component.
163 Example:
164 \begin{python}
165 x=Symbol('x', (2,2))
166 y=3*x
167 print y.coeff(x)
168 print y.coeff(x[1,1])
169 \end{python}
170 will print:
171 \begin{python}
172 [[3 3]
173 [3 3]]
174 [[0 0]
175 [0 3]]
176 \end{python}
177 \end{methoddesc}
178 \begin{methoddesc}[Symbol]{diff}{symbols}
179 Takes the derivative of the \SYMBOL object which the method belongs to with respect to the symbols specified in the argument symbols.
180 \end{methoddesc}
181 \begin{methoddesc}[Symbol]{evalf}{}
182 Applies the sympy.evalf operation on all elements of the \SYMBOL which are of type or inherit from sympy.Basic.
183 \end{methoddesc}
184 \begin{methoddesc}[Symbol]{expand}{}
185 Applies the sympy.expand operation on all elements in this \SYMBOL.
186 \end{methoddesc}
187 %\begin{methoddesc}[Symbol]{getDataSubstitutions}{}
188 %end{methoddesc}
189 \begin{methoddesc}[Symbol]{getDim}{}
190 Returns the symbol's spatial dimensionality, or -1 if undefined.
191 \end{methoddesc}
192 \begin{methoddesc}[Symbol]{getRank}{}
193 Returns the symbol's rank which is equal to the length of the shape.
194 \end{methoddesc}
195 \begin{methoddesc}[Symbol]{getShape}{}
196 Returns the shape of this \SYMBOL.
197 \end{methoddesc}
198 \pagebreak
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}
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}
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 \pagebreak
283 \begin{methoddesc}[NonlinearPDE]{setOptions}{opts}
284 \begin{verbatim}
285 Allows setting options for the nonlinear PDE.
286 The supported options are:
287 tolerance
288 error tolerance for the Newton method
289 iteration_steps_max
290 maximum number of Newton iterations
291 omega_min
292 minimum relaxation factor
293 atol
294 solution norms less than atol are assumed to be atol.
295 This can be useful if one of your solutions is expected to
296 be zero.
297 quadratic_convergence_limit
298 if the norm of the Newton-Raphson correction is reduced by
299 less than quadratic_convergence_limit between two iteration
300 steps, quadratic convergence is assumed.
301 simplified_newton_limit
302 if the norm of the defect is reduced by less than
303 simplified_newton_limit between two iteration steps and
304 quadratic convergence is detected, the iteration switches to the
305 simplified Newton-Raphson scheme.
306 \end{verbatim}
308 \end{methoddesc}
309 \begin{methoddesc}[NonlinearPDE]{setValue}{
310 \optional{X}\optional{, Y}
311 \optional{, y}\optional{, y_contact}
312 \optional{, y_dirac}
313 \optional{, q}\optional{, r}}
314 Assigns new values to coefficients. By default all values are assumed to be
315 zero\footnote{In fact, it is assumed they are not present by assigning the
316 value \code{escript.Data()}. This can be used by the solver library to reduce
317 computational costs.}.
318 If the new coefficient value is not a \Data object, it is converted into a
319 \Data object in the appropriate \FunctionSpace.
320 \end{methoddesc}
322 \subsection{Symconsts class}
323 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.
324 \begin{verbatim}
325 usage:
326 symconsts.pi this provides a Symbol object
328 Available constants currently include:
329 pi and e
330 \end{verbatim}

  ViewVC Help
Powered by ViewVC 1.1.26