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

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}