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

revision 2974 by artak, Wed Mar 3 04:27:00 2010 UTC revision 2975 by ahallam, Thu Mar 4 05:44:12 2010 UTC
# Line 30  We use \file{onedheatdiff001b.py} as the Line 30  We use \file{onedheatdiff001b.py} as the
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
# Line 49  model = Rectangle(l0=mx,l1=my,n0=ndx, n1 Line 47  model = Rectangle(l0=mx,l1=my,n0=ndx, n1
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}
# Line 72  We can easily apply the backward Euler s Line 69  We can easily apply the backward Euler s
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^{(n-1)}  \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)}
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
# Line 93  The fundamental difference to the 1D cas Line 90  The fundamental difference to the 1D cas
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 \verb|kronecker| returns this matrix;  \esc function \verb|kronecker| returns this matrix;
96  \begin{equation}  \begin{equation}
# Line 125  this means $\frac{\partial T^{(n)}}{\par Line 122 this means$\frac{\partial T^{(n)}}{\par
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,
# Line 155  kappa = kappai * whereNegative(bound) + Line 152  kappa = kappai * whereNegative(bound) +
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 \verb|Function| function space as expected for the  Notice that we are using the sample points of the \verb|Function| function space as expected for the
155  PDE coefficient \verb|A|\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.}  PDE coefficient \verb|A|\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()
# Line 168  T= Ti*whereNegative(bound)+Tc*(1-whereNe Line 165  T= Ti*whereNegative(bound)+Tc*(1-whereNe
165  where \verb|Ti| and \verb|Tc| are the initial temperature  where \verb|Ti| and \verb|Tc| 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 \verb|x| and \verb|bound| to refer to the appropriate  notice that we reset \verb|x| and \verb|bound| 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}}
# Line 183  sample points of a PDE solution\footnote Line 180  sample points of a PDE solution\footnote
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
# Line 202  function in other scripts more easily. Line 199  function in other scripts more easily.
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 \verb|linspace|  .  are the dimensions in the $x$ and $y$ directions) which is done using the \modnumpy function \verb|linspace|  .
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
# Line 213  yi = np.linspace(0.0,my,75) Line 210  yi = np.linspace(0.0,my,75)
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.
# Line 245  of packages available. Here we use the p Line 242  of packages available. Here we use the p
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

Legend:
 Removed from v.2974 changed lines Added in v.2975