/[escript]/trunk/doc/cookbook/example03.tex
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2952 by gross, Thu Feb 25 09:52:10 2010 UTC revision 2969 by artak, Wed Mar 3 04:27:00 2010 UTC
# Line 20  Line 20 
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 cross-section as shown in \reffig{fig:twodhdmodel}. We will implement the same  we model a cross-section 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
# Line 31  We use \file{onedheatdiff001b.py} as the Line 31  We use \file{onedheatdiff001b.py} as the
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
# Line 49  model = Rectangle(l0=mx,l1=my,n0=ndx, n1 Line 49  model = Rectangle(l0=mx,l1=my,n0=ndx, n1
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
# Line 89  $f = q \hackscore{H} + \frac{\rho c\hack Line 89  $f = q \hackscore{H} + \frac{\rho c\hack
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
# Line 111  mypde.setValue(A=kappa*kronecker(model), Line 111  mypde.setValue(A=kappa*kronecker(model),
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}
# Line 126  where we have  $(n\hackscore{0},n\hacksc Line 126  where we have  $(n\hackscore{0},n\hacksc
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 sub-domain.  curved interface between the two sub-domain.
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   (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}   (x-x\hackscore{0})^{2}+(y-y\hackscore{0})^{2} = s^{2}
# Line 147  Defining these conditions within the scr Line 147  Defining these conditions within the scr
147  \end{python}  \end{python}
148  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| command again to help define the material constants and temperatures in our domain.
149  If \verb|kappai| and \verb|kappac| are the  If \verb|kappai| and \verb|kappac| 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(x-ic)-r  bound = length(x-ic)-r
# Line 167  T= Ti*whereNegative(bound)+Tc*(1-whereNe Line 167  T= Ti*whereNegative(bound)+Tc*(1-whereNe
167  \end{python}  \end{python}
168  where \verb|Ti| and \verb|Tc| are the initial temperature  where \verb|Ti| and \verb|Tc| 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 \verb|x| and \verb|bound| to refer to the appropriate  notice that we reset \verb|x| and \verb|bound| 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]
# Line 213  yi = np.linspace(0.0,my,75) Line 213  yi = np.linspace(0.0,my,75)
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)
# Line 239  The results for selected time steps are Line 239  The results for selected time steps are
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).
# Line 280  created by \verb|saveVTK|. Line 280  created by \verb|saveVTK|.
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}

Legend:
Removed from v.2952  
changed lines
  Added in v.2969

  ViewVC Help
Powered by ViewVC 1.1.26