# Diff of /trunk/doc/cookbook/example03.tex

revision 2658 by ahallam, Thu Sep 10 02:58:44 2009 UTC revision 2948 by gross, Thu Feb 25 04:54:30 2010 UTC
# Line 1  Line 1
1
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2009 by University of Queensland  % Copyright (c) 2003-2010 by University of Queensland
5  % Earth Systems Science Computational Center (ESSCC)  % Earth Systems Science Computational Center (ESSCC)
6  % http://www.uq.edu.au/esscc  % http://www.uq.edu.au/esscc
7  %  %
# Line 11  Line 11
11  %  %
12  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14  \section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion}  \begin{figure}[t]
\sslist{twodheatdiff001.py and cblib.py}
%\label{Sec:2DHD}
Building upon our success from the 1D models it is now prudent to expand our domain by another dimension. For this example we will be using a very simple magmatic intrusion as the basis for our model. The simulation will be a single event where some molten granite has formed a hemisphericle dome at the base of some cold sandstone country rock. A hemisphere is symmetric so taking a cross-section through its centre will effectively model a 3D problem in 2D. New concepts will include non-linear boundaries and the ability to prescribe location specific variables.

\begin{figure}[h!]
15  \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}  \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}
16  \caption{2D model: granitic intrusion of sandstone country rock.}  \caption{2D model: granitic intrusion of sandstone country rock.}
17  \label{fig:twodhdmodel}  \label{fig:twodhdmodel}
18  \end{figure}  \end{figure}
19
20  To expand upon our 1D problem, the domain must first be expanded. This will be done in our definition phase by creating a square domain in $x$ and $y$ that is 600 meters along each side \reffig{fig:twodhdmodel}. The number of discrete spatial cells will be 100. The radius of the intrusion will be 200 meters  And the location of the centre of the intrusion will be at the 300 meter mark on the x-axis. The domain variables are;  \sslist{twodheatdiff001.py and cblib.py}
21  \begin{verbatim}
22    Building upon our success from the 1D models, it is now prudent to expand our domain by another dimension.
23    For this example we will be using a very simple magmatic intrusion as the basis for our model. The simulation will be a single event where some molten granite has formed a cylindrical dome at the base of some cold sandstone country rock. Assuming that the cylinder is very long
24    we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement the same
25    diffusion model as we have use for the granite blocks in \refSec{Sec:1DHDv00}
26    but will add the second spatial dimension and show how to define
27    variables depending on the location of the domain.
28    We use \file{onedheatdiff001b.py} as the starting point for develop this model.
29
30    \section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion}
31    \label{Sec:2DHD}
32
33    To expand upon our 1D problem, the domain must first be expanded. In fact, we have solved a two dimensional problem already but didn't put much
34    attention to the second dimension. This will be changed now.
35    In our definition phase by creating a square domain in $x$ and $y$\footnote{In \esc the notation
36    $x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. The number of discrete spatial cells will be 100 in either direction. The radius of the intrusion will be $200$ meters with the centre located at the $300$ meter mark on the $x$-axis. The domain variables are;
37    \begin{python}
38  mx = 600*m #meters - model length  mx = 600*m #meters - model length
39  my = 600*m #meters - model width  my = 600*m #meters - model width
40  ndx = 100 #mesh steps in x direction  ndx = 150 #mesh steps in x direction
41  ndy = 100 #mesh steps in y direction  ndy = 150 #mesh steps in y direction
42  r = 200*m #meters - radius of intrusion  r = 200*m #meters - radius of intrusion
43  ic = [300, 0] #centre of intrusion (meters)  ic = [300*m, 0] #coordinates of the centre of intrusion (meters)
44  q=0.*Celsius #our heat source temperature is now zero  qH=0.*J/(sec*m**3) #our heat source temperature is zero
45  \end{verbatim}  \end{python}
46  The next step is to define our variables for each material in the model in a manner similar to the previous tutorial. Note that each material has its own unique set of values. The time steps and set up for the domain remain as before. Prior to setting up the PDE the boundary between the two materials must be established. The distance $s$ between two points in car Cartesian coordinates is defined as:  As before we use
47    \begin{python}
48    model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
49    \end{python}
50    to generate the domain.
51
52    There are two fundamental changes to the PDE has we have discussed PDEs in \refSec{Sec:1DHDv00}. Firstly,
53    as the material coefficients for granite and sandstone is different, we need to deal with
54    PDE coefficients which vary with there location in the domain. Secondly, we need
55    to deal with the second spatial dimension. We will look at these two modification at the same time.
56    In fact our temperature model \refEq{eqn:hd} we have used in \refSec{Sec:1DHDv00} applies for the
57    1D case with constant material parameter only. For the more general case we are interested
58    in this chapter the correct model equation is
59
60   (x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2} = s^{2}  \rho c\hackscore p \frac{\partial T}{\partial t} -  \frac{\partial }{\partial x} \kappa \frac{\partial T}{\partial x} -  \frac{\partial }{\partial y} \kappa \frac{\partial T}{\partial y} = q\hackscore H
61    \label{eqn:hd2}
62
63  If $[x_{0},y_{0}]$ is the point $c$ the centre of the semi-circle that defines our intrusion then for all the points $[x,y]$ in our solution space we can define a distance to $c$. Hence, and points that fall within the radius $r$ of our intrusion will have a corresponding value $s < r$ and all those belonging to the country rock will have a value $s > r$. By subtracting $r$ from both of these conditions we find $s-r < 0$ for all intrusion points and $s-r > 0$ for all country rock points. Defining these conditions within the script is quite simple and is done using the following command:  Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and
64  \begin{verbatim}  for the case of constant material parameters $\frac{\partial }{\partial x} \kappa = \kappa \frac{\partial }{\partial x}$ so this new equation coincides with simplified model equation for this case. It is more convenient
65    to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla};
66    \label{eqn:Tform nabla}
67    \rho c\hackscore p \frac{\partial T}{\partial t}
68    -\nabla \cdot \kappa \nabla T = q\hackscore H
69
70    We can easily apply the backward Euler scheme as before to end up with
71
72    \frac{\rho c\hackscore p}{h} T^{(n)} -\nabla \cdot \kappa \nabla T^{(n)}  = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}
73    \label{eqn:hdgenf2}
74
75    which is very similar to \refEq{eqn:hdgenf} used to define the temperature update in the simple 1D case.
76    The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under
77    the light of the more general case we need to revisit the \esc PDE form as
78    shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients
79    the form needs to be read as
80    \label{eqn:commonform2D 2}
81    -\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x}
82    -\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y}
83    -\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x}
84    -\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y}
85    + Du = f
86
87    So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and
88    $f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}$ as we have used before (see \refEq{ESCRIPT SET}) we need to set
89    \label{eqn: kappa general}
90    A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0
91
92    The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$
93    which brings in the second dimension. Important to notice that the fact that the coefficients
94    of the PDE may depend on their location in the domain now does not influence the usage of the PDE form. This is very convenient as we can introduce spatial dependence to the PDE coefficients without modification to the way we call the PDE solver.
95
96    A very convenient way to define the matrix $A$ is required in \refEq{eqn: kappa general} is using the
97    Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The
98    \esc function \verb|kronecker| returns this matrix;
99
100    \verb|kronecker(model)| = \left[
101    \begin{array}{cc}
102     1 & 0 \\
103     0 & 1 \\
104    \end{array}
105    \right]
106
107    As the argument \verb|model| represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix
108    (In case of a three-dimensional domain a $3 \times 3$ matrix is returned). The call
109    \begin{python}
110    mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
111    \end{python}
112    sets the PDE coefficients according to \refEq{eqn: kappa general}.
113
114    Before we turn the question how we set $\kappa$ we need to check the boundary conditions. As
115    pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case
116    this assumptions translates to;
117
118    -n \cdot \kappa \nabla T^{(n)} =
119    -n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} - n\hackscore{1} \cdot  \kappa \frac{\partial T^{(n)}}{\partial y} = 0
120    \label{eq:hom flux}
121
122    which sets the normal component of the heat flux $- \kappa \cdot (\frac{\partial T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means that the regions is insulated which is what we want.
123    On the left and right face of the domain where we have $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$
124    this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain
125    where we have  $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$.
126
127    \section{Setting Variable PDE Coefficients}
128    Now we need to look into the problem how we define the material coefficients
129    $\kappa$ and $\rho c\hackscore p$ depending on there location in the domain.
130    We have used the technique we discuss here already when we set up the initial
131    temperature in the granite block example in \refSec{Sec:1DHDv00}. However,
132    the situation is more complicated here as we have to deal with a
133    curved interface between the two sub-domain.
134
135    Prior to setting up the PDE the interface between the two materials must be established.
136    The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as:
137
138     (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}
139
140    If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes the centre of the semi-circle of our intrusion, then for all the points $[x,y]$ in our model we can calculate a distance to $ic$.
141    All the points that fall within the given radius $r$ of our intrusion will have a corresponding
142    value $s < r$ and all those belonging to the country rock will have a value $s > r$. By subtracting $r$ from both of these conditions we find $s-r < 0$ for all intrusion points and $s-r > 0$
143    for all country rock points.
144    Defining these conditions within the script is quite simple and is done using the following command:
145    \begin{python}
146   bound = length(x-ic)-r #where the boundary will be located   bound = length(x-ic)-r #where the boundary will be located
147  \end{verbatim}  \end{python}
148  This definition of the boundary can now be used with the \verb wherePositive()  and \verb whereNegative()  commands from before to help define the material constants and temperatures in our domain. By examining the general form we solved in the earlier tutorials, it is obvious that both \verb A  and \verb D  depend on the predefined variables. To set these variables accordingly and complete our PDE we use:  This definition of the boundary can now be used with \verb|whereNegative| command again to help define the material constants and temperatures in our domain.
149  \begin{verbatim}  If \verb|kappai| and \verb|kappac| are the
150  A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound)  thermal conductivities for the intrusion material granite and for the surrounding sandstone we set;
151  D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound)  \begin{python}
152    x=Function(model).getX()
153    bound = length(x-ic)-r
154    kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound))
155    mypde.setValue(A=kappa*kronecker(model))
156    \end{python}
157    Notice that we are using the sample points of the \verb|Function| function space as expected for the
158    PDE coefficient \verb|A|\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.}
159    The corresponding statements are used to set $\rho c\hackscore p$.
160
mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc)
\end{verbatim}
161  Our PDE has now been properly established. The initial conditions for temperature are set out in a similar matter:  Our PDE has now been properly established. The initial conditions for temperature are set out in a similar matter:
162  \begin{verbatim}  \begin{python}
163  #defining the initial temperatures.  #defining the initial temperatures.
164   T= Ti*whereNegative(bound)+Tc*wherePositive(bound)  x=Solution(model).getX()
165  \end{verbatim}  bound = length(x-ic)-r
166  The iteration process now begins as before, but using our new conditions for \verb D  as defined above.  T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
167    \end{python}
168  \subsection{Contouring escript data}  where \verb|Ti| and \verb|Tc| are the initial temperature
169  It is possible to contour our solution using \modmpl . Unfortunately the \modmpl contouring function only accepts regularly gridded data. As our solution is not regularly gridded, it is necessary to interpolate our solution onto a regular grid. First we extract the model coordinates using \verb getX  these are then transformed to a numpy array and into individual $x$ and $y$ arrays. We also need to generate our regular grid which is done using the \modnumpy function \verb linspace  .  in the regions of the granite and surrounding sandstone, respectively. It is important to
170  \begin{verbatim}  notice that we have reset \verb|x| and \verb|bound| to refer to the appropriate
171  # rearrage mymesh to suit solution function space for contouring        sample points of a PDE solution\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}.
172  oldspacecoords=model.getX()
173  coords=Data(oldspacecoords, T.getFunctionSpace())  \begin{figure}[h]
174  coordX, coordY = toXYTuple(coords)  \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}}
175    \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}}
176    \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}}
177    \caption{2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.}
178    \label{fig:twodhdans}
179    \end{figure}
180
181    \section{Contouring \esc data using \modmpl.}
182    \label{Sec:2DHD plot}
183    To complete our transition from a 1D to a 2D model we also need to modify the
184    plotting procedure. As before we use the  \modmpl to do the plotting
185    in this case a contour plot. For 2D contour plots \modmpl expects that the
186    data are regularly gridded. We have no control on the location and ordering of the sample points
187    used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid.
188
189    In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as
190    \verb|numpy| array to hand the values to \modmpl. This can easily be extended to extract both the
191    $x$ and the $y$ coordinates;
192    \begin{python}
193    import numpy as np
194    def toXYTuple(coords):
195        coords = np.array(coords.toListOfTuples()) #convert to Tuple
196        coordX = coords[:,0] #X components.
197        coordY = coords[:,1] #Y components.
198        return coordX,coordY
199    \end{python}
200    For convenience we have put this function into \file{clib.py} file so we can use this
201    function in other scripts more easily.
202
203
204    We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$
205    are the dimensions in $x$ and $y$ direction) which is done using the \modnumpy function \verb|linspace|  .
206    \begin{python}
207    from clib import toXYTuple
208    # get sample points for temperature as  for contouring
209    coordX, coordY = toXYTuple(T.getFunctionSpace().getX())
210  # create regular grid  # create regular grid
211  xi = np.linspace(0.0,mx,100)  xi = np.linspace(0.0,mx,75)
212  yi = np.linspace(0.0,my,100)  yi = np.linspace(0.0,my,75)
213  \end{verbatim}  \end{python}
214  The remainder of our contouring commands reside within the \verb while  loop so that a new contour is generated for each time step. For each time step the solution much be regridded for \modmpl using the \verb griddata  function. This function interprets an irregular grid and solution from \verb tempT  , \verb xi   and \verb yi  this is transformed to the new coordinates defined by \verb coordX  and \verb coordY  with an output \verb zi  . It is now possible to use the \verb contourf  function which generates a colour filled contour. The colour gradient of our plots is set via the command \verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page. Our results are then contoured, visually adjusted using the \modmpl functions and then saved to file. \verb pl.clf()  clears the figure in readiness for the next time iteration.  The values \verb|[xi[k], yi[l]]| are the grid points.
215
216  \begin{verbatim}  The remainder of our contouring commands reside within a \verb while  loop so that a new contour is generated for each time step. For each time step the solution must be regridded for \modmpl using the \verb griddata  function. This function interprets a potentially irregularly located values \verb tempT  at locations defined by \verb coordX  and \verb coordY as values at the new coordinates of a rectangular grid defined by
217    \verb xi   and \verb yi . The output is \verb zi  . It is now possible to use the \verb contourf  function which generates colour filled contours. The colour gradient of our plots is set via the command \verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page\footnote{see \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured, visually adjusted using the \modmpl functions and then saved to file. \verb pl.clf()  clears the figure in readiness for the next time iteration.
218    \begin{python}
219  #grid the data.  #grid the data.
220  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
221  # contour the gridded data, plotting dots at the randomly spaced data points.  # contour the gridded data, plotting dots at the randomly spaced data points.
# Line 81  pl.xlabel("Horizontal Displacement (m)") Line 229  pl.xlabel("Horizontal Displacement (m)")
229  pl.ylabel("Depth (m)")  pl.ylabel("Depth (m)")
230  pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)  pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)
231  pl.clf()          pl.clf()
232  \end{verbatim}  \end{python}
233    The function \verb|pl.contour| is used to add labeled contour lines to the plot.
234    The results for selected time steps are shown in \reffig{fig:twodhdans}.
235
\begin{figure}[h!]
\centerline{\includegraphics[width=4.in]{figures/heatrefraction050}}
\caption{2D model: Total temperature distribution ($T$) at time $t=50$.}
\label{fig:twodhdmodel}
\end{figure}
236
237
239
240    \sslist{twodheatdiffvtk.py}
241    An alternative approach to \modmpl for visualization is the usage of a package which base on
242    visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There is a variety
243    of package available. Here we will use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example.
244
245    \mayavi is an interactive, GUI driven tool which is
246    really designed to visualize large three dimensional data sets where \modmpl
247    is not suitable. But it is very useful when it comes to two dimensional problems.
248    The decision which tool is best is finally the user's decision. The main
249    difference between using \mayavi (and other VTK based tools)
250    or \modmpl is the fact that actually visualization is detached from the
251    calculation by writing the results to external files
252    and import them into \mayavi. In 3D where the best camera position for rendering a scene is not obvious
253    before the results are available. Therefore the user may need to try
254    different position before the best is found. Moreover, in many cases in 3D the interactive
255    visualization is the only way to really understand the results (e.g. using stereographic projection).
256
257    To write the temperatures at each time step to data files in the VTK file format one
258    needs to insert a \verb|saveVTK| call into the code;
259    \begin{python}
260    while t<=tend:
261          i+=1 #counter
262          t+=h #current time
263          mypde.setValue(Y=qH+T*rhocp/h)
264          T=mypde.getSolution()
265          saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T)
266    \end{python}
267    The data files, eg. \file{data.001.vtu}, contains all necessary information to
268    visualize the temperature and can directly processed by \mayavi. Notice that there is no
269    regridding required. It is recommended to use the file extension \file{.vtu} for files
270    created by \verb|saveVTK|.
271
272    \begin{figure}[ht]
273    \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}}
274    \caption{\mayavi start up Window.}
275    \label{fig:mayavi window}
276    \end{figure}
277
278    \begin{figure}[ht]
279    \centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}}
280    \caption{\mayavi data control window.}
281    \label{fig:mayavi window2}
282    \end{figure}
283    After you have run the script you will find the
284    result files \file{data.*.vtu} in the result directory \file{data/twodheatdiff}. Run the
285    command
286    \begin{python}
287    >> mayavi2 -d data.001.vtu -m Surface &
288    \end{python}
289    from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}.
290    The right hand side shows the temperature at the first time step. To show
291    the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)}
292    at the top left hand side. This will bring up a new window similar to tye window shown in  \reffig{fig:mayavi window2}. By clicking at the arrows in the top right corner you can move forwards and backwards in time.
293    You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The
294    name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used
295    in the \verb|saveVTK| call. In this menu item you can select other results
296    you may have written to the output file for visualization.
297
298    \textbf{For the advanced user}: Using the \modmpl to visualize spatially distributed data
299    is not MPI compatible. However, the \verb|saveVTK| function can be used with MPI. In fact,
300    the function collects the values of the sample points spread across processor ranks into a single.
301    For more details on writing scripts for parallel computing please consult the \emph{user's guide}.

Legend:
 Removed from v.2658 changed lines Added in v.2948