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

revision 2975 by ahallam, Thu Mar 4 05:44:12 2010 UTC revision 2979 by ahallam, Tue Mar 9 02:54:32 2010 UTC
# Line 19  Line 19
19
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
23  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  domain by another dimension.
24  we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement the same  For this example we use a very simple magmatic intrusion as the basis for our
25    model. The simulation will be a single event where some molten granite has
26    formed a cylindrical dome at the base of some cold sandstone country rock.
27    Assuming that the cylinder is very long
28    we model a cross-section as shown in \reffig{fig:twodhdmodel}. We will implement
29    the same
30  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}
31  but will add the second spatial dimension and show how to define  but will add the second spatial dimension and show how to define
32  variables depending on the location of the domain.  variables depending on the location of the domain.
33  We use \file{onedheatdiff001b.py} as the starting point for develop this model.  We use \file{onedheatdiff001b.py} as the starting point for develop this model.
34
35  \section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic Intrusion}  \section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic
36    Intrusion}
37  \label{Sec:2DHD}  \label{Sec:2DHD}
38
39  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  To expand upon our 1D problem, the domain must first be expanded. In fact, we
40  $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;  have solved a two dimensional problem already but essentially ignored the second
41    dimension. In our definition phase we create a square domain in $x$ and
42    $y$\footnote{In \esc the notation
43    $x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.}
44    that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the
45    number of discrete spatial cells to 150 in both direction and the radius of the
46    intrusion to $200$ meters with the centre located at the $300$ meter mark on the
47    $x$-axis. Thus, the domain variables are;
48  \begin{python}  \begin{python}
49  mx = 600*m #meters - model length  mx = 600*m #meters - model length
50  my = 600*m #meters - model width  my = 600*m #meters - model width
# Line 47  model = Rectangle(l0=mx,l1=my,n0=ndx, n1 Line 60  model = Rectangle(l0=mx,l1=my,n0=ndx, n1
60  \end{python}  \end{python}
61  to generate the domain.  to generate the domain.
62
63  There are two fundamental changes to the PDE that we have discussed in \refSec{Sec:1DHDv00}. Firstly,  There are two fundamental changes to the PDE that we have discussed in
64  because the material coefficients for granite and sandstone are different, we need to deal with  \refSec{Sec:1DHDv00}. Firstly,
65  PDE coefficients which vary with their location in the domain. Secondly, we need  because the material coefficients for granite and sandstone are different, we
66  to deal with the second spatial dimension. We can investigate these two modification at the same time.  need to deal with
67  In fact, the temperature model \refEq{eqn:hd} we utilised in \refSec{Sec:1DHDv00} applied for the  PDE coefficients which vary with their location in the domain. Secondly, we
68  1D case with a constant material parameter only. For the more general case examined in this chapter, the correct model equation is  need
69    to deal with the second spatial dimension. We can investigate these two
70    modification at the same time.
71    In fact, the temperature model \refEq{eqn:hd} we utilised in
72    \refSec{Sec:1DHDv00} applied for the
73    1D case with a constant material parameter only. For the more general case
74    examined in this chapter, the correct model equation is
75
76  \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}
77    \kappa \frac{\partial T}{\partial x} -  \frac{\partial }{\partial y} \kappa
78    \frac{\partial T}{\partial y} = q\hackscore H
79  \label{eqn:hd2}  \label{eqn:hd2}
80
81  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
82  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  for the case of constant material parameters $\frac{\partial }{\partial x} 83 to write this equation using the$\nabla$notation as we have already seen in \refEq{eqn:commonform nabla}; \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides
84    with a simplified model equation for this case. It is more convenient
85    to write this equation using the $\nabla$ notation as we have already seen in
86    \refEq{eqn:commonform nabla};
87  \label{eqn:Tform nabla}  \label{eqn:Tform nabla}
88  \rho c\hackscore p \frac{\partial T}{\partial t}  \rho c\hackscore p \frac{\partial T}{\partial t}
89  -\nabla \cdot \kappa \nabla T = q\hackscore H  -\nabla \cdot \kappa \nabla T = q\hackscore H
90
91  We can easily apply the backward Euler scheme as before to end up with  We can easily apply the backward Euler scheme as before to end up with
92
93  \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)}  =
94    q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}
95  \label{eqn:hdgenf2}  \label{eqn:hdgenf2}
96
97  which is very similar to \refEq{eqn:hdgenf} used to define the temperature solution in the simple 1D case.  which is very similar to \refEq{eqn:hdgenf} used to define the temperature
98  The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under  update in the simple 1D case.
99    The difference is in the second order derivate term $\nabla \cdot \kappa \nabla 100 T^{(n)}$. Under
101  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
102  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
103  the form needs to be read as  the form needs to be read as
# Line 82  the form needs to be read as Line 109  the form needs to be read as
109  + Du = f  + Du = f
110
111  So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and  So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and
112  $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  $f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}$ as we have used
113    before (see \refEq{ESCRIPT SET}) we need to set
114  \label{eqn: kappa general}  \label{eqn: kappa general}
115  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
116
117  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
118    to zero but $\kappa$,
119  which brings in the second dimension. Important to notice that the coefficients  which brings in the second dimension. Important to notice that the coefficients
120  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
121    usage of the PDE form. This is very convenient as we can introduce spatial
122  A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can be carried out using the  dependence to the PDE coefficients without modification to the way we call the
123  Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The  PDE solver.
124
125    A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can
126    be carried out using the
127    Kronecker $\delta$ symbol\footnote{see
128    \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The
129  \esc function \verb|kronecker| returns this matrix;  \esc function \verb|kronecker| returns this matrix;
130
131  \verb|kronecker(model)| = \left[  \verb|kronecker(model)| = \left[
# Line 101  Kronecker $\delta$ symbol\footnote{see \ Line 135  Kronecker $\delta$ symbol\footnote{see \
135  \end{array}  \end{array}
136  \right]  \right]
137
138  As the argument \verb|model| represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix  As the argument \verb|model| represents a two dimensional domain the matrix is
139  (In case of a three-dimensional domain a $3 \times 3$ matrix is returned). The call  returned as $2 \times 2$ matrix
140    (In case of a three-dimensional domain a $3 \times 3$ matrix is returned). The
141    call
142  \begin{python}  \begin{python}
143  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)  mypde.setValue(A=kappa*kronecker(model),D=rhocp/h)
144  \end{python}  \end{python}
145  sets the PDE coefficients according to \refEq{eqn: kappa general}.    sets the PDE coefficients according to \refEq{eqn: kappa general}.
146
147  We need to check the boundary conditions before we turn to the question: how we set $\kappa$. As  We need to check the boundary conditions before we turn to the question: how we
148  pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case  set $\kappa$. As
149    pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary
150    conditions. In our case
151  this assumptions translates to;  this assumptions translates to;
152
153  -n \cdot \kappa \nabla T^{(n)} =  -n \cdot \kappa \nabla T^{(n)} =
154  -n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} - n\hackscore{1} \cdot  \kappa \frac{\partial T^{(n)}}{\partial y} = 0  -n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x} -
155    n\hackscore{1} \cdot  \kappa \frac{\partial T^{(n)}}{\partial y} = 0
156  \label{eq:hom flux}  \label{eq:hom flux}
157
158  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.  which sets the normal component of the heat flux $- \kappa \cdot (\frac{\partial 159 On the left and right face of the domain where we have$(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means
160  this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain  that the regions is insulated which is what we want.
161  where we have  $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$.  On the left and right face of the domain where we have
162    $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$
163    this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on
164    the domain
165    where we have  $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is
166    $\frac{\partial T^{(n)}}{\partial y}=0$.
167
168  \section{Setting Variable PDE Coefficients}  \section{Setting Variable PDE Coefficients}
169  Now we need to look into the problem of how we define the material coefficients  Now we need to look into the problem of how we define the material coefficients
170  $\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.
171  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
172    \refSec{Sec:1DHDv00}
173  to set up the initial temperature. However,  to set up the initial temperature. However,
174  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
175  curved interface between the two sub-domain.  curved interface between the two sub-domain.
176
177  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
178  The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as:  established.
179    The distance $s\ge 0$ between two points $[x,y]$ and
180    $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as:
181
182   (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}   (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}
183
184  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$.  If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes
185  All the points that fall within the given radius $r$ of our intrusion will have a corresponding  the centre of the semi-circle of our intrusion, then for all the points $[x,y]$
186  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$  in our model we can calculate a distance to $ic$.
187    All the points that fall within the given radius $r$ of our intrusion will have
188    a corresponding
189    value $s < r$ and all those belonging to the country rock will have a value $s > 190 r$. By subtracting $r$ from both of these conditions we find $s-r < 0$ for all
191    intrusion points and $s-r > 0$
192  for all country rock points.  for all country rock points.
193  Defining these conditions within the script is quite simple and is done using the following command:  Defining these conditions within the script is quite simple and is done using
194    the following command:
195  \begin{python}  \begin{python}
196   bound = length(x-ic)-r #where the boundary will be located   bound = length(x-ic)-r #where the boundary will be located
197  \end{python}  \end{python}
198  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.  This definition of the boundary can now be used with \verb|whereNegative|
199    command again to help define the material constants and temperatures in our
200    domain.
201  If \verb|kappai| and \verb|kappac| are the  If \verb|kappai| and \verb|kappac| are the
202  thermal conductivities for the intrusion material granite and for the surrounding sandstone, then we set;  thermal conductivities for the intrusion material granite and for the
203    surrounding sandstone, then we set;
204  \begin{python}  \begin{python}
205  x=Function(model).getX()  x=Function(model).getX()
206  bound = length(x-ic)-r  bound = length(x-ic)-r
207  kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound))  kappa = kappai * whereNegative(bound) + kappac * (1-whereNegative(bound))
208  mypde.setValue(A=kappa*kronecker(model))  mypde.setValue(A=kappa*kronecker(model))
209  \end{python}  \end{python}
210  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
211  PDE coefficient \verb|A|\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.}  as expected for the
212    PDE coefficient \verb|A|\footnote{For the experienced user: use
213    \texttt{x=mypde.getFunctionSpace("A").getX()}.}
214  The corresponding statements are used to set $\rho c\hackscore p$.  The corresponding statements are used to set $\rho c\hackscore p$.
215
216  Our PDE has now been properly established. The initial conditions for temperature are set out in a similar manner:  Our PDE has now been properly established. The initial conditions for
217    temperature are set out in a similar manner:
218  \begin{python}  \begin{python}
219  #defining the initial temperatures.  #defining the initial temperatures.
220  x=Solution(model).getX()  x=Solution(model).getX()
# Line 163  bound = length(x-ic)-r Line 222  bound = length(x-ic)-r
222  T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))  T= Ti*whereNegative(bound)+Tc*(1-whereNegative(bound))
223  \end{python}  \end{python}
224  where \verb|Ti| and \verb|Tc| are the initial temperature  where \verb|Ti| and \verb|Tc| are the initial temperature
225  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
226    important to
227  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
228  sample points of a PDE solution\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}.  sample points of a PDE solution\footnote{For the experienced user: use
229    \texttt{x=mypde.getFunctionSpace("r").getX()}.}.
230
231  \begin{figure}[ht]  \begin{figure}[ht]
232  \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}}  \centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}}
233  \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}}  \centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}}
234  \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}}  \centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}}
235  \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.}  \caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step
236    $1$, $20$ and $200$. Contour lines show temperature.}
237  \label{fig:twodhdans}  \label{fig:twodhdans}
238  \end{figure}  \end{figure}
239
# Line 180  sample points of a PDE solution\footnote Line 242  sample points of a PDE solution\footnote
242  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
243  plotting procedure. As before we use the  \modmpl to do the plotting  plotting procedure. As before we use the  \modmpl to do the plotting
244  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
245  data are regularly gridded. We have no control over the location and ordering of the sample points  data are regularly gridded. We have no control over the location and ordering of
246  used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid.  the sample points
247    used to represent the solution. Therefore it is necessary to interpolate our
248  In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as  solution onto a regular grid.
249  \verb|numpy| array to hand the values to \modmpl. This can easily be extended to extract both the
250    In \refSec{sec: plot T} we have already learned how to extract the $x$
251    coordinates of sample points as
252    \verb|numpy| array to hand the values to \modmpl. This can easily be extended to
253    extract both the
254  $x$ and the $y$ coordinates;  $x$ and the $y$ coordinates;
255  \begin{python}  \begin{python}
256  import numpy as np  import numpy as np
# Line 194  def toXYTuple(coords): Line 260  def toXYTuple(coords):
260      coordY = coords[:,1] #Y components.      coordY = coords[:,1] #Y components.
261      return coordX,coordY      return coordX,coordY
262  \end{python}  \end{python}
263  For convenience we have put this function into \file{clib.py} file so we can use this  For convenience we have put this function into \file{clib.py} file so we can use
264    this
265  function in other scripts more easily.  function in other scripts more easily.
266
267
268  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$
269  are the dimensions in the $x$ and $y$ directions) which is done using the \modnumpy function \verb|linspace|  .  are the dimensions in the $x$ and $y$ directions) which is done using the
270    \modnumpy function \verb|linspace|  .
271  \begin{python}  \begin{python}
272  from clib import toXYTuple  from clib import toXYTuple
273  # get sample points for temperature as  for contouring        # get sample points for temperature as  for contouring
# Line 210  yi = np.linspace(0.0,my,75) Line 278  yi = np.linspace(0.0,my,75)
278  \end{python}  \end{python}
279  The values \verb|[xi[k], yi[l]]| are the grid points.  The values \verb|[xi[k], yi[l]]| are the grid points.
280
281  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  The remainder of our contouring commands reside within a \verb while  loop so
282  \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.  that a new contour is generated for each time step. For each time step the
283    solution must be regridded for \modmpl using the \verb griddata  function. This
284    function interprets irregularly located values \verb tempT  at locations defined
285    by \verb coordX  and \verb coordY  as values at the new coordinates of a
286    rectangular grid defined by
287    \verb xi   and \verb yi . The output is \verb zi . It is now possible to use the
288    \verb contourf  function which generates colour filled contours. The colour
289    gradient of our plots is set via the command
290    \verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page\footnote{see
291    \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured,
292    visually adjusted using the \modmpl functions and then saved to file. \verb
293    pl.clf()  clears the figure in readiness for the next time iteration.
294  \begin{python}  \begin{python}
295  #grid the data.  #grid the data.
296  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)  zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi)
# Line 227  pl.ylabel("Depth (m)") Line 306  pl.ylabel("Depth (m)")
306  pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i)  pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i)
307  pl.clf()          pl.clf()
308  \end{python}  \end{python}
309  The function \verb|pl.contour| is used to add labeled contour lines to the plot.  The function \verb|pl.contour| is used to add labeled contour lines to the
310    plot.
311  The results for selected time steps are shown in \reffig{fig:twodhdans}.  The results for selected time steps are shown in \reffig{fig:twodhdans}.
312
313
# Line 235  The results for selected time steps are Line 315  The results for selected time steps are
316
317  \sslist{example03b.py}  \sslist{example03b.py}
318  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
319  visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There are variety  which base on
320  of packages available. Here we use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example.  visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}.
321    There are variety
322    of packages available. Here we use the package \mayavi\footnote{see
323    \url{http://code.enthought.com/projects/mayavi/}} as an example.
324
325  \mayavi is an interactive, GUI driven tool which is  \mayavi is an interactive, GUI driven tool which is
326  really designed to visualize large three dimensional data sets where \modmpl  really designed to visualize large three dimensional data sets where \modmpl
327  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
328  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)  problems.
329    The decision of which tool is the best can be subjective and the user should
330    determine which package they require and are most comfortable with. The main
331    difference between using \mayavi (and other VTK based tools)
332  or \modmpl is that the actual visualization is detached from the  or \modmpl is that the actual visualization is detached from the
333  calculation by writing the results to external files  calculation by writing the results to external files
334  and importing 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
335    scene is not obvious
336  before the results are available. Therefore the user may need to try  before the results are available. Therefore the user may need to try
337  different position before the best is found. Moreover, in many cases a 3D interactive  different position before the best is found. Moreover, in many cases a 3D
338  visualization is the only way to really understand the results (e.g. using stereographic projection).  interactive
339    visualization is the only way to really understand the results (e.g. using
340    stereographic projection).
341
342  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
343    one
344  needs to insert a \verb|saveVTK| call into the code;  needs to insert a \verb|saveVTK| call into the code;
345  \begin{python}  \begin{python}
346  while t<=tend:  while t<=tend:
# Line 261  while t<=tend: Line 351  while t<=tend:
351        saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T)        saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T)
352  \end{python}  \end{python}
353  The data files, eg. \file{data.001.vtu}, contains all necessary information to  The data files, eg. \file{data.001.vtu}, contains all necessary information to
354  visualize the temperature and can directly processed by \mayavi. Notice that there is no  visualize the temperature and can directly processed by \mayavi. Notice that
355  regridding required. It is recommended to use the file extension \file{.vtu} for files  there is no
356    regridding required. It is recommended to use the file extension \file{.vtu} for
357    files
358  created by \verb|saveVTK|.  created by \verb|saveVTK|.
359
360  \begin{figure}[ht]  \begin{figure}[ht]
# Line 277  created by \verb|saveVTK|. Line 369  created by \verb|saveVTK|.
369  \label{fig:mayavi window2}  \label{fig:mayavi window2}
370  \end{figure}  \end{figure}
371  After you run the script you will find the  After you run the script you will find the
372  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}.
373    Run the
374  command  command
375  \begin{python}  \begin{python}
376  >> mayavi2 -d data.001.vtu -m Surface &  >> mayavi2 -d data.001.vtu -m Surface &
377  \end{python}  \end{python}
378  from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}.  from the result directory. \mayavi will start up a window similar to
379    \reffig{fig:mayavi window}.
380  The right hand side shows the temperature at the first time step. To show  The right hand side shows the temperature at the first time step. To show
381  the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)}  the results at other time steps you can click at the item \texttt{VTK XML file
382  at the top left hand side. This will bring up a new window similar to the 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.  (data.001.vtu) (timeseries)}
383  You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The  at the top left hand side. This will bring up a new window similar to the window
384  name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used  shown in  \reffig{fig:mayavi window2}. By clicking at the arrows in the top
385    right corner you can move forwards and backwards in time.
386    You will also notice the text \textbf{T} next to the item \texttt{Point scalars
387    name}. The
388    name \textbf{T} corresponds to the keyword argument name \texttt{T} we have
389    used
390  in the \verb|saveVTK| call. In this menu item you can select other results  in the \verb|saveVTK| call. In this menu item you can select other results
391  you may have written to the output file for visualization.  you may have written to the output file for visualization.
392
\textbf{For the advanced user}: Using the \modmpl to visualize spatially distributed data
is not MPI compatible. However, the \verb|saveVTK| function can be used with MPI. In fact,
the function collects the values of the sample points spread across processor ranks into a single.
For more details on writing scripts for parallel computing please consult the \emph{user's guide}.
393    \textbf{For the advanced user}: Using the \modmpl to visualize spatially
394    distributed data
395    is not MPI compatible. However, the \verb|saveVTK| function can be used with
396    MPI. In fact,
397    the function collects the values of the sample points spread across processor
398    ranks into a single.
399    For more details on writing scripts for parallel computing please consult the
400    \emph{user's guide}.

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