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

revision 2632 by ahallam, Wed Aug 26 22:18:19 2009 UTC revision 2645 by ahallam, Thu Sep 3 02:20:33 2009 UTC
# Line 13  Line 13
13
14  \section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion}  \section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion}
15  %\label{Sec:2DHD}  %\label{Sec:2DHD}
16   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 semetric 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 abilitity to prescribe location specific variables.   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.
17
18  \begin{figure}[h!]  \begin{figure}[h!]
19  \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}  \centerline{\includegraphics[width=4.in]{figures/twodheatdiff}}
# Line 21  Line 21
21  \label{fig:twodhdmodel}  \label{fig:twodhdmodel}
22  \end{figure}  \end{figure}
23
24  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 eacy 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;  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;
25  \begin{verbatim}  \begin{verbatim}
26  mx = 600 # model lenght  mx = 600*m #meters - model length
27  my = 600 # model width  my = 600*m #meters - model width
28  ndx = 100 # steps in x direction  ndx = 100 #mesh steps in x direction
29  ndy = 100 # steps in y direction  ndy = 100 #mesh steps in y direction
30  r = 200 # radius of intrusion  r = 200*m #meters - radius of intrusion
31  ic = [300, 0] #centre of intrusion  ic = [300, 0] #centre of intrusion (meters)
32    q=0.*Celsius #our heat source temperature is now zero
33  \end{verbatim}  \end{verbatim}
34  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:  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:
35
36   (x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2} = s^{2}   (x_{1}-x_{0})^{2}+(y_{1}-y_{0})^{2} = s^{2}
37
38  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 coresponding 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 into the code is quite simple and is done using the following command:  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:
39  \begin{verbatim}  \begin{verbatim}
40   bound = length(x-ic)-r #where the boundary will be located   bound = length(x-ic)-r #where the boundary will be located
41  \end{verbatim}  \end{verbatim}
# Line 49  Our PDE has now been properly establishe Line 50  Our PDE has now been properly establishe
50  \begin{verbatim}  \begin{verbatim}
51   T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures.   T= Ti*whereNegative(bound)+Tc*wherePositive(bound) #defining the initial temperatures.
52  \end{verbatim}  \end{verbatim}
53  The itteration process now begins as before, but using our new conditions for \verb D  as defined above.  The iteration process now begins as before, but using our new conditions for \verb D  as defined above.
54
55  \subsection{Contouring escript data}  \subsection{Contouring escript data}
56    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  .
57  \begin{verbatim}  \begin{verbatim}
58  # rearrage mymesh to suit solution function space        # rearrage mymesh to suit solution function space for contouring
59  oldspacecoords=model.getX()  oldspacecoords=model.getX()
60  coords=Data(oldspacecoords, T.getFunctionSpace())  coords=Data(oldspacecoords, T.getFunctionSpace())
61  coords = np.array(coords.toListOfTuples())  coordX, coordY = toXYTuple(coords)
coordX = coords[:,0]
coordY = coords[:,1]
62  # create regular grid  # create regular grid
63  xi = np.linspace(0.0,600.0,100)  xi = np.linspace(0.0,mx,100)
64  yi = np.linspace(0.0,600.0,100)  yi = np.linspace(0.0,my,100)
65  \end{verbatim}  \end{verbatim}
66    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.
67
68    \begin{verbatim}
69    #grid the data.
70    zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
71    # contour the gridded data, plotting dots at the randomly spaced data points.
72    pl.matplotlib.pyplot.autumn()
73    pl.contourf(xi,yi,zi,10)
74    CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
75    pl.clabel(CS, inline=1, fontsize=8)
76    pl.axis([0,600,0,600])
77    pl.title("Heat diffusion from an intrusion.")
78    pl.xlabel("Horizontal Displacement (m)")
79    pl.ylabel("Depth (m)")
80    pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i)
81    pl.clf()
82    \end{verbatim}
83
84    \begin{figure}[h!]
85    \centerline{\includegraphics[width=4.in]{figures/heatrefraction050}}
86    \caption{2D model: Total temperature distribution ($T$) at time $t=50$.}
87    \label{fig:twodhdmodel}
88    \end{figure}

Legend:
 Removed from v.2632 changed lines Added in v.2645