/[escript]/trunk/doc/cookbook/example03.tex
ViewVC logotype

Contents of /trunk/doc/cookbook/example03.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 17577 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2018 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Apache License, version 2.0
8 % http://www.apache.org/licenses/LICENSE-2.0
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development 2012-2013 by School of Earth Sciences
12 % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 %
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
16 \begin{figure}[t]
17 \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}
18 \caption{Example 3: 2D model: granitic intrusion of sandstone country rock}
19 \label{fig:twodhdmodel}
20 \end{figure}
21
22 \sslist{example03a.py and cblib.py}
23
24 Building upon our success from the 1D models, it is now prudent to expand our
25 domain by another dimension.
26 For this example we use a very simple magmatic intrusion as the basis for our
27 model. The simulation will be a single event where some molten granite has
28 formed a cylindrical dome at the base of some cold sandstone country rock.
29 Assuming that the cylinder is very long
30 we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement
31 the same
32 diffusion model as we have used for the granite blocks in \refSec{Sec:1DHDv00}
33 but will add the second spatial dimension and show how to define
34 variables depending on the location of the domain.
35 We use \file{onedheatdiff001b.py} as the starting point to develop this model.
36
37 \section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic
38 Intrusion}
39 \label{Sec:2DHD}
40
41 To expand upon our 1D problem, the domain must first be expanded. In fact, we
42 have solved a two dimensional problem already but essentially ignored the second
43 dimension. In our definition phase we create a square domain in $x$ and
44 $y$\footnote{In \esc the notation
45 $x_{0}$ and $x_{1}$ is used for $x$ and $y$, respectively.}
46 that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the
47 number of discrete spatial cells to $150$ in both directions and the radius of
48 the intrusion to $200$ meters with the centre located at the $300$ meter mark
49 on the $x$-axis. Thus, the domain variables are;
50 \begin{python}
51 mx = 600*m #meters - model length
52 my = 600*m #meters - model width
53 ndx = 150 #mesh steps in x direction
54 ndy = 150 #mesh steps in y direction
55 r = 200*m #meters - radius of intrusion
56 ic = [300*m, 0] #coordinates of the centre of intrusion (meters)
57 qH=0.*J/(sec*m**3) #our heat source temperature is zero
58 \end{python}
59 As before we use
60 \begin{python}
61 model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
62 \end{python}
63 to generate the domain.
64
65 There are two fundamental changes to the PDE that we have discussed in
66 \refSec{Sec:1DHDv00}. Firstly,
67 because the material coefficients for granite and sandstone are different, we
68 need to deal with
69 PDE coefficients which vary with their location in the domain. Secondly, we
70 need to deal with the second spatial dimension. We can investigate these two
71 modifications at the same time.
72 In fact, the temperature model \refEq{eqn:hd} we utilised in
73 \refSec{Sec:1DHDv00} applied for the
74 1D case with a constant material parameter only. For the more general case
75 examined in this chapter, the correct model equation is
76 \begin{equation}
77 \rho c_p \frac{\partial T}{\partial t} - \frac{\partial }{\partial x}
78 \kappa \frac{\partial T}{\partial x} - \frac{\partial }{\partial y} \kappa
79 \frac{\partial T}{\partial y} = q_H
80 \label{eqn:hd2}
81 \end{equation}
82 Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and
83 for the case of constant material parameters $\frac{\partial }{\partial x}
84 \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides
85 with a simplified model equation for this case. It is more convenient
86 to write this equation using the $\nabla$ notation as we have already seen in
87 \refEq{eqn:commonform nabla};
88 \begin{equation}\label{eqn:Tform nabla}
89 \rho c_p \frac{\partial T}{\partial t}
90 -\nabla \cdot \kappa \nabla T = q_H
91 \end{equation}
92 We can easily apply the backward Euler scheme as before to end up with
93 \begin{equation}
94 \frac{\rho c_p}{h} T^{(n)} -\nabla \cdot \kappa \nabla T^{(n)} =
95 q_H + \frac{\rho c_p}{h} T^{(n-1)}
96 \label{eqn:hdgenf2}
97 \end{equation}
98 which is very similar to \refEq{eqn:hdgenf} used to define the temperature
99 update in the simple 1D case.
100 The difference is in the second order derivative term
101 $\nabla \cdot \kappa \nabla T^{(n)}$. Under the light of the more general case
102 we need to revisit the \esc PDE form as shown in \ref{eqn:commonform2D}.
103 For the 2D case with variable PDE coefficients the form needs to be read as
104 \begin{equation}\label{eqn:commonform2D 2}
105 -\frac{\partial }{\partial x} A_{00}\frac{\partial u}{\partial x}
106 -\frac{\partial }{\partial x} A_{01}\frac{\partial u}{\partial y}
107 -\frac{\partial }{\partial y} A_{10}\frac{\partial u}{\partial x}
108 -\frac{\partial }{\partial x} A_{11}\frac{\partial u}{\partial y}
109 + Du = f
110 \end{equation}
111 So besides the settings $u=T^{(n)}$, $D = \frac{\rho c _{p}}{h}$ and
112 $f = q _{H} + \frac{\rho c_p}{h} T^{(n-1)}$ as we have used
113 before (see \refEq{ESCRIPT SET}) we need to set
114 \begin{equation}\label{eqn: kappa general}
115 A_{00}=A_{11}=\kappa; A_{01}=A_{10}=0
116 \end{equation}
117 The fundamental difference to the 1D case is that $A_{11}$ is not set
118 to zero but $\kappa$,
119 which brings in the second dimension. It is important to note that the
120 coefficients of the PDE may depend on their location in the domain which does
121 not influence the usage of the PDE form. This is very convenient as we can
122 introduce spatial dependence to the PDE coefficients without modification to
123 the way we call the PDE solver.
124
125 A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general}
126 can be carried out using the Kronecker $\delta$ symbol\footnote{see
127 \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The
128 \esc function \verb|kronecker| returns this matrix;
129 \begin{equation}
130 \verb|kronecker(model)| = \left[
131 \begin{array}{cc}
132 1 & 0 \\
133 0 & 1 \\
134 \end{array}
135 \right]
136 \end{equation}
137 As the argument \verb|model| represents a two dimensional domain the matrix is
138 returned as a $2 \times 2$ matrix
139 (in the case of a three-dimensional domain a $3 \times 3$ matrix is returned).
140 The call
141 \begin{python}
142 mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
143 \end{python}
144 sets the PDE coefficients according to \refEq{eqn: kappa general}.
145
146 We need to check the boundary conditions before we turn to the question: how do
147 we set $\kappa$. As pointed out \refEq{NEUMAN 2} makes certain assumptions on
148 the boundary conditions. In our case these assumptions translate to;
149 \begin{equation}
150 -n \cdot \kappa \nabla T^{(n)} =
151 -n_{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} -
152 n_{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0
153 \label{eq:hom flux}
154 \end{equation}
155 which sets the normal component of the heat flux $- \kappa \cdot (\frac{\partial
156 T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means
157 that the region is insulated which is what we want.
158 On the left and right face of the domain where we have
159 $(n_{0},n_{1} ) = (\mp 1,0)$
160 this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on
161 the domain
162 where we have $(n_{0},n_{1} ) = (\pm 1,0)$ this is
163 $\frac{\partial T^{(n)}}{\partial y}=0$.
164
165 \section{Setting variable PDE Coefficients}
166 Now we need to look into the problem of how we define the material coefficients
167 $\kappa$ and $\rho c_p$ depending on their location in the domain.
168 We can make use of the technique used in the granite block example in
169 \refSec{Sec:1DHDv00}
170 to set up the initial temperature. However,
171 the situation is more complicated here as we have to deal with a
172 curved interface between the two sub-domains.
173
174 Prior to setting up the PDE, the interface between the two materials must be
175 established.
176 The distance $s\ge 0$ between two points $[x,y]$ and
177 $[x_{0},y_{0}]$ in Cartesian coordinates is defined as:
178 \begin{equation}
179 (x-x_{0})^{2}+(y-y_{0})^{2} = s^{2}
180 \end{equation}
181 If we define the point $[x_{0},y_{0}]$ as $ic$ which denotes
182 the centre of the semi-circle of our intrusion, then for all the points $[x,y]$
183 in our model we can calculate a distance to $ic$.
184 All the points that fall within the given radius $r$ of our intrusion will have
185 a corresponding
186 value $s < r$ and all those belonging to the country rock will have a value $s >
187 r$. By subtracting $r$ from both of these conditions we find $s-r < 0$ for all
188 intrusion points and $s-r > 0$ for all country rock points.
189 Defining these conditions within the script is quite simple and is done using
190 the following command:
191 \begin{python}
192 bound = length(x-ic)-r #where the boundary will be located
193 \end{python}
194 This definition of the boundary can now be used with the \verb|whereNegative|
195 command again to help define the material constants and temperatures in our
196 domain.
197 If \verb|kappai| and \verb|kappac| are the
198 thermal conductivities for the intrusion material granite and for the
199 surrounding sandstone, then we set;
200 \begin{python}
201 x=Function(model).getX()
202 bound = length(x-ic)-r
203 kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound))
204 mypde.setValue(A=kappa*kronecker(model))
205 \end{python}
206 Notice that we are using the sample points of the \verb|Function| function space
207 as expected for the PDE coefficient \verb|A|\footnote{For the experienced user: use
208 \texttt{x=mypde.getFunctionSpace("A").getX()}.}.
209 The corresponding statements are used to set $\rho c_p$.
210
211 Our PDE has now been properly established. The initial conditions for
212 temperature are set out in a similar manner:
213 \begin{python}
214 #defining the initial temperatures.
215 x=Solution(model).getX()
216 bound = length(x-ic)-r
217 T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
218 \end{python}
219 where \verb|Ti| and \verb|Tc| are the initial temperature in the regions of the
220 granite and surrounding sandstone, respectively. It is important to
221 notice that we reset \verb|x| and \verb|bound| to refer to the appropriate
222 sample points of a PDE solution\footnote{For the experienced user: use
223 \texttt{x=mypde.getFunctionSpace("r").getX()}.}.
224
225 \begin{figure}[ht]
226 \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}}
227 \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}}
228 \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}}
229 \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step
230 $1$, $20$ and $200$. Contour lines show temperature.}
231 \label{fig:twodhdans}
232 \end{figure}
233
234 \section{Contouring \esc Data using \modmpl.}
235 \label{Sec:2DHD plot}
236 To complete our transition from a 1D to a 2D model we also need to modify the
237 plotting procedure. As before we use \modmpl to do the plotting in this case a
238 contour plot. For 2D contour plots \modmpl expects that the data are regularly
239 gridded. We have no control over the location and ordering of the sample points
240 used to represent the solution. Therefore it is necessary to interpolate our
241 solution onto a regular grid.
242
243 In \refSec{sec: plot T} we have already learned how to extract the $x$
244 coordinates of sample points as
245 \verb|numpy| array to hand the values to \modmpl. This can easily be extended
246 to extract both the $x$ and the $y$ coordinates;
247 \begin{python}
248 import numpy as np
249 def toXYTuple(coords):
250 coords = np.array(coords.toListOfTuples()) #convert to Tuple
251 coordX = coords[:,0] #X components.
252 coordY = coords[:,1] #Y components.
253 return coordX,coordY
254 \end{python}
255 For convenience we have put this function into \file{clib.py} so we can use
256 this function in other scripts more easily.
257
258 We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$
259 are the dimensions in the $x$ and $y$ directions) which is done using the
260 \modnumpy function \verb|linspace|.
261 \begin{python}
262 from clib import toXYTuple
263 # get sample points for temperature as for contouring
264 coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
265 # create regular grid
266 xi = np.linspace(0.0,mx,75)
267 yi = np.linspace(0.0,my,75)
268 \end{python}
269 The values \verb|[xi[k], yi[l]]| are the grid points.
270
271 The remainder of our contouring commands resides within a \verb|while| loop so
272 that a new contour is generated for each time step. For each time step the
273 solution must be re-gridded for \modmpl using the \verb|griddata| function. This
274 function interprets irregularly located values \verb|tempT| at locations defined
275 by \verb|coordX| and \verb|coordY| as values at the new coordinates of a
276 rectangular grid defined by
277 \verb|xi| and \verb|yi|. The output is \verb|zi|. It is now possible to use the
278 \verb|contourf| function which generates colour filled contours. The colour
279 gradient of our plots is set via the command
280 \verb|pl.matplotlib.pyplot.autumn()|, other colours are listed on the \modmpl web page\footnote{see
281 \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured,
282 visually adjusted using the \modmpl functions and then saved to a file.
283 \verb|pl.clf()| clears the figure in readiness for the next time iteration.
284 \begin{python}
285 #grid the data.
286 zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
287 # contour the gridded data, plotting dots at the randomly spaced data points.
288 pl.matplotlib.pyplot.autumn()
289 pl.contourf(xi,yi,zi,10)
290 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
291 pl.clabel(CS, inline=1, fontsize=8)
292 pl.axis([0,600,0,600])
293 pl.title("Heat diffusion from an intrusion.")
294 pl.xlabel("Horizontal Displacement (m)")
295 pl.ylabel("Depth (m)")
296 pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i)
297 pl.clf()
298 \end{python}
299 The function \verb|pl.contour| is used to add labelled contour lines to the
300 plot.
301 The results for selected time steps are shown in \reffig{fig:twodhdans}.
302
303
304 \section{Advanced Visualisation using VTK}
305
306 \sslist{example03b.py}
307 An alternative approach to \modmpl for visualisation is the usage of a package
308 which is based on the Visualization Toolkit (VTK) library\footnote{see \url{http://www.vtk.org/}}.
309 There is a variety of packages available. Here we use the package \mayavi\footnote{see
310 \url{http://code.enthought.com/projects/mayavi/}} as an example.
311
312 \mayavi is an interactive, GUI driven tool which is
313 really designed to visualise large three dimensional data sets where \modmpl
314 is not suitable. But it is also very useful when it comes to two dimensional
315 problems.
316 The decision of which tool is the best can be subjective and users should
317 determine which package they require and are most comfortable with. The main
318 difference between using \mayavi (or other VTK based tools) and \modmpl is that
319 the actual visualisation is detached from the calculation by writing the
320 results to external files and importing them into \mayavi. In 3D the best
321 camera position for rendering a scene is not obvious before the results are
322 available. Therefore the user may need to try different settings before the
323 best is found. Moreover, in many cases a 3D interactive visualisation is the
324 only way to really understand the results (e.g. using stereographic projection).
325
326 To write the temperatures at each time step to data files in the VTK file format
327 one needs to import \verb|saveVTK| from the \weipa module and call it:
328 \begin{python}
329 from esys.weipa import saveVTK
330 while t<=tend:
331 i+=1 #counter
332 t+=h #current time
333 mypde.setValue(Y=qH+T*rhocp/h)
334 T=mypde.getSolution()
335 saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T)
336 \end{python}
337 The data files, e.g. \file{data.001.vtu}, contain all necessary information to
338 visualise the temperature and can be directly processed by \mayavi. Note that
339 there is no re-gridding required. The file extension \file{.vtu} is
340 automatically added if not supplied to \verb|saveVTK|.
341
342 \begin{figure}[ht]
343 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}}
344 \caption{Example 3b: \mayavi start up window}
345 \label{fig:mayavi window}
346 \end{figure}
347
348 \begin{figure}[ht]
349 \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}}
350 \caption{Example 3b: \mayavi data control window}
351 \label{fig:mayavi window2}
352 \end{figure}
353 After you run the script you will find the
354 result files \file{data.*.vtu} in the result directory \file{data/example03}.
355 Run the command
356 \begin{python}
357 >> mayavi2 -d data.001.vtu -m Surface &
358 \end{python}
359 from the result directory. \mayavi will start up a window similar to
360 \reffig{fig:mayavi window}.
361 The right hand side shows the temperature at the first time step. To show
362 the results at other time steps you can click at the item \texttt{VTK XML file
363 (data.001.vtu) (timeseries)}
364 at the top left hand side. This will bring up a new window similar to the window
365 shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top
366 right corner you can move forwards and backwards in time.
367 You will also notice the text \textbf{T} next to the item \texttt{Point scalars
368 name}. The
369 name \textbf{T} corresponds to the keyword argument name \texttt{T} we have
370 used in the \verb|saveVTK| call. In this menu item you can select other results
371 you may have written to the output file for visualisation.
372
373 \textbf{For the advanced user}: Using the \modmpl to visualise spatially
374 distributed data is not MPI compatible. However, the \verb|saveVTK| function
375 can be used with MPI. In fact, the function collects the values of the sample
376 points spread across processor ranks into a single file.
377 For more details on writing scripts for parallel computing please consult the
378 \emph{user's guide}.
379

  ViewVC Help
Powered by ViewVC 1.1.26