30 
\section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic Intrusion} 
\section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic Intrusion} 
31 
\label{Sec:2DHD} 
\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 
To expand upon our 1D problem, the domain must first be expanded. In fact, we have solved a two dimensional problem already but essentially ignored the second dimension. In our definition phase we create a square domain in $x$ and $y$\footnote{In \esc the notation 

attention to the second dimension. 


In our definition phase we create a square domain in $x$ and $y$\footnote{In \esc the notation 

34 
$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the number of discrete spatial cells to 150 in both direction and the radius of the intrusion to $200$ meters with the centre located at the $300$ meter mark on the $x$axis. Thus, the domain variables are; 
$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the number of discrete spatial cells to 150 in both direction and the radius of the intrusion to $200$ meters with the centre located at the $300$ meter mark on the $x$axis. Thus, the domain variables are; 
35 
\begin{python} 
\begin{python} 
36 
mx = 600*m #meters  model length 
mx = 600*m #meters  model length 
47 
\end{python} 
\end{python} 
48 
to generate the domain. 
to generate the domain. 
49 


50 
There are two fundamental changes to the PDE that we have discussed PDEs in \refSec{Sec:1DHDv00}. Firstly, 
There are two fundamental changes to the PDE that we have discussed in \refSec{Sec:1DHDv00}. Firstly, 
51 
as the material coefficients for granite and sandstone is different, we need to deal with 
because the material coefficients for granite and sandstone are different, we need to deal with 
52 
PDE coefficients which vary with their location in the domain. Secondly, we need 
PDE coefficients which vary with their location in the domain. Secondly, we need 
53 
to deal with the second spatial dimension. We look at these two modification at the same time. 
to deal with the second spatial dimension. We can investigate these two modification at the same time. 
54 
In fact our temperature model \refEq{eqn:hd} we have used in \refSec{Sec:1DHDv00} applies for the 
In fact, the temperature model \refEq{eqn:hd} we utilised in \refSec{Sec:1DHDv00} applied for the 
55 
1D case with constant material parameter only. For the more general case we are interested 
1D case with a constant material parameter only. For the more general case examined in this chapter, the correct model equation is 

in this chapter the correct model equation is 

56 
\begin{equation} 
\begin{equation} 
57 
\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 
\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 
58 
\label{eqn:hd2} 
\label{eqn:hd2} 
59 
\end{equation} 
\end{equation} 
60 
Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
61 
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 
for the case of constant material parameters $\frac{\partial }{\partial x} \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides with a simplified model equation for this case. It is more convenient 
62 
to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 
to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 
63 
\begin{equation}\label{eqn:Tform nabla} 
\begin{equation}\label{eqn:Tform nabla} 
64 
\rho c\hackscore p \frac{\partial T}{\partial t} 
\rho c\hackscore p \frac{\partial T}{\partial t} 
69 
\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
70 
\label{eqn:hdgenf2} 
\label{eqn:hdgenf2} 
71 
\end{equation} 
\end{equation} 
72 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature update in the simple 1D case. 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature solution in the simple 1D case. 
73 
The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 
The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 
74 
the light of the more general case we need to revisit the \esc PDE form as 
the light of the more general case we need to revisit the \esc PDE form as 
75 
shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
90 
which brings in the second dimension. Important to notice that the coefficients 
which brings in the second dimension. Important to notice that the coefficients 
91 
of the PDE may depend on their location in the domain and 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. 
of the PDE may depend on their location in the domain and 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. 
92 


93 
A very convenient way to define the matrix $A$ is required in \refEq{eqn: kappa general} is using the 
A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can be carried out using the 
94 
Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
95 
\esc function \verbkronecker returns this matrix; 
\esc function \verbkronecker returns this matrix; 
96 
\begin{equation} 
\begin{equation} 
122 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 
123 


124 
\section{Setting Variable PDE Coefficients} 
\section{Setting Variable PDE Coefficients} 
125 
Now we need to look into the problem how we define the material coefficients 
Now we need to look into the problem of how we define the material coefficients 
126 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
127 
We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 
We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 
128 
to set up the initial temperature. However, 
to set up the initial temperature. However, 
152 
mypde.setValue(A=kappa*kronecker(model)) 
mypde.setValue(A=kappa*kronecker(model)) 
153 
\end{python} 
\end{python} 
154 
Notice that we are using the sample points of the \verbFunction function space as expected for the 
Notice that we are using the sample points of the \verbFunction function space as expected for the 
155 
PDE coefficient \verbA\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 
PDE coefficient \verbA\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 
156 
The corresponding statements are used to set $\rho c\hackscore p$. 
The corresponding statements are used to set $\rho c\hackscore p$. 
157 


158 
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 manner: 
159 
\begin{python} 
\begin{python} 
160 
#defining the initial temperatures. 
#defining the initial temperatures. 
161 
x=Solution(model).getX() 
x=Solution(model).getX() 
165 
where \verbTi and \verbTc are the initial temperature 
where \verbTi and \verbTc are the initial temperature 
166 
in the regions of the granite and surrounding sandstone, respectively. It is important to 
in the regions of the granite and surrounding sandstone, respectively. It is important to 
167 
notice that we reset \verbx and \verbbound to refer to the appropriate 
notice that we reset \verbx and \verbbound to refer to the appropriate 
168 
sample points of a PDE solution\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
sample points of a PDE solution\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
169 


170 
\begin{figure}[ht] 
\begin{figure}[ht] 
171 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
180 
To complete our transition from a 1D to a 2D model we also need to modify the 
To complete our transition from a 1D to a 2D model we also need to modify the 
181 
plotting procedure. As before we use the \modmpl to do the plotting 
plotting procedure. As before we use the \modmpl to do the plotting 
182 
in this case a contour plot. For 2D contour plots \modmpl expects that the 
in this case a contour plot. For 2D contour plots \modmpl expects that the 
183 
data are regularly gridded. We have no control on the location and ordering of the sample points 
data are regularly gridded. We have no control over the location and ordering of the sample points 
184 
used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid. 
used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid. 
185 


186 
In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as 
In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as 
199 


200 


201 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
202 
are the dimensions in $x$ and $y$ direction) which is done using the \modnumpy function \verblinspace . 
are the dimensions in the $x$ and $y$ directions) which is done using the \modnumpy function \verblinspace . 
203 
\begin{python} 
\begin{python} 
204 
from clib import toXYTuple 
from clib import toXYTuple 
205 
# get sample points for temperature as for contouring 
# get sample points for temperature as for contouring 
210 
\end{python} 
\end{python} 
211 
The values \verb[xi[k], yi[l]] are the grid points. 
The values \verb[xi[k], yi[l]] are the grid points. 
212 


213 
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 
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 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 
214 
\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. 
\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. 
215 
\begin{python} 
\begin{python} 
216 
#grid the data. 
#grid the data. 
242 
\mayavi is an interactive, GUI driven tool which is 
\mayavi is an interactive, GUI driven tool which is 
243 
really designed to visualize large three dimensional data sets where \modmpl 
really designed to visualize large three dimensional data sets where \modmpl 
244 
is not suitable. But it is very useful when it comes to two dimensional problems. 
is not suitable. But it is very useful when it comes to two dimensional problems. 
245 
The decision which tool is the best depends on the user's decision. The main 
The decision of which tool is the best can be subjective and the user should determine which package they require and are most comfortable with. The main difference between using \mayavi (and other VTK based tools) 
246 
difference between using \mayavi (and other VTK based tools) 
or \modmpl is that the actual visualization is detached from the 

or \modmpl is the fact that actually visualization is detached from the 

247 
calculation by writing the results to external files 
calculation by writing the results to external files 
248 
and import them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 
and importing them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 
249 
before the results are available. Therefore the user may need to try 
before the results are available. Therefore the user may need to try 
250 
different position before the best is found. Moreover, in many cases in 3D the interactive 
different position before the best is found. Moreover, in many cases a 3D interactive 
251 
visualization is the only way to really understand the results (e.g. using stereographic projection). 
visualization is the only way to really understand the results (e.g. using stereographic projection). 
252 


253 
To write the temperatures at each time step to data files in the VTK file format one 
To write the temperatures at each time step to data files in the VTK file format one 