20 
\sslist{example03a.py and cblib.py} 
\sslist{example03a.py and cblib.py} 
21 


22 
Building upon our success from the 1D models, it is now prudent to expand our domain by another dimension. 
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 
For this example we use 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 crosssection as shown in \reffig{fig:twodhdmodel}. We will implement the same 
we model a crosssection 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} 
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 
but will add the second spatial dimension and show how to define 
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 didn't put much 
34 
attention to the second dimension. This will be changed now. 
attention to the second dimension. 
35 
In our definition phase by creating a square domain in $x$ and $y$\footnote{In \esc the notation 
In our definition phase we create 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; 
$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; 
37 
\begin{python} 
\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 
49 
\end{python} 
\end{python} 
50 
to generate the domain. 
to generate the domain. 
51 


52 
There are two fundamental changes to the PDE has we have discussed PDEs in \refSec{Sec:1DHDv00}. Firstly, 
There are two fundamental changes to the PDE that 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 
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 
PDE coefficients which vary with their 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. 
to deal with the second spatial dimension. We 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 
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 
1D case with constant material parameter only. For the more general case we are interested 
58 
in this chapter the correct model equation is 
in this chapter the correct model equation is 
89 
\begin{equation}\label{eqn: kappa general} 
\begin{equation}\label{eqn: kappa general} 
90 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
91 
\end{equation} 
\end{equation} 
92 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$ 
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 
which brings in the second dimension. Important to notice 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. 
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. 
95 


96 
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$ 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 
Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
111 
\end{python} 
\end{python} 
112 
sets the PDE coefficients according to \refEq{eqn: kappa general}. 
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 
We need to check the boundary conditions before we turn to the question: how we set $\kappa$. As 
115 
pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 
pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 
116 
this assumptions translates to; 
this assumptions translates to; 
117 
\begin{equation} 
\begin{equation} 
126 


127 
\section{Setting Variable PDE Coefficients} 
\section{Setting Variable PDE Coefficients} 
128 
Now we need to look into the problem how we define the material coefficients 
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. 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
130 
We have used the technique we discuss here already when we set up the initial 
We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 
131 
temperature in the granite block example in \refSec{Sec:1DHDv00}. However, 
to set up the initial temperature. However, 
132 
the situation is more complicated here as we have to deal with a 
the situation is more complicated here as we have to deal with a 
133 
curved interface between the two subdomain. 
curved interface between the two subdomain. 
134 


135 
Prior to setting up the PDE the interface between the two materials must be established. 
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: 
The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
137 
\begin{equation} 
\begin{equation} 
138 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
147 
\end{python} 
\end{python} 
148 
This definition of the boundary can now be used with \verbwhereNegative command again to help define the material constants and temperatures in our domain. 
This definition of the boundary can now be used with \verbwhereNegative command again to help define the material constants and temperatures in our domain. 
149 
If \verbkappai and \verbkappac are the 
If \verbkappai and \verbkappac are the 
150 
thermal conductivities for the intrusion material granite and for the surrounding sandstone we set; 
thermal conductivities for the intrusion material granite and for the surrounding sandstone, then we set; 
151 
\begin{python} 
\begin{python} 
152 
x=Function(model).getX() 
x=Function(model).getX() 
153 
bound = length(xic)r 
bound = length(xic)r 
167 
\end{python} 
\end{python} 
168 
where \verbTi and \verbTc are the initial temperature 
where \verbTi and \verbTc are the initial temperature 
169 
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 
170 
notice that we have reset \verbx and \verbbound to refer to the appropriate 
notice that we reset \verbx and \verbbound to refer to the appropriate 
171 
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 experience user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
172 


173 
\begin{figure}[ht] 
\begin{figure}[ht] 
213 
\end{python} 
\end{python} 
214 
The values \verb[xi[k], yi[l]] are the grid points. 
The values \verb[xi[k], yi[l]] are the grid points. 
215 


216 
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 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. 
\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} 
\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) 
239 


240 
\sslist{example03b.py} 
\sslist{example03b.py} 
241 
An alternative approach to \modmpl for visualization is the usage of a package which base on 
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 
visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There are variety 
243 
of package available. Here we will use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example. 
of packages available. Here we 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 
\mayavi is an interactive, GUI driven tool which is 
246 
really designed to visualize large three dimensional data sets where \modmpl 
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. 
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 
The decision which tool is the best depends on the user's decision. The main 
249 
difference between using \mayavi (and other VTK based tools) 
difference between using \mayavi (and other VTK based tools) 
250 
or \modmpl is the fact that actually visualization is detached from the 
or \modmpl is the fact that actually visualization is detached from the 
251 
calculation by writing the results to external files 
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 
and import them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 
253 
before the results are available. Therefore the user may need to try 
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 
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). 
visualization is the only way to really understand the results (e.g. using stereographic projection). 
280 
\caption{Example 3b: \mayavi data control window.} 
\caption{Example 3b: \mayavi data control window.} 
281 
\label{fig:mayavi window2} 
\label{fig:mayavi window2} 
282 
\end{figure} 
\end{figure} 
283 
After you have run the script you will find the 
After you run the script you will find the 
284 
result files \file{data.*.vtu} in the result directory \file{data/example03}. Run the 
result files \file{data.*.vtu} in the result directory \file{data/example03}. Run the 
285 
command 
command 
286 
\begin{python} 
\begin{python} 