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

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

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

revision 2971 by gross, Thu Feb 25 09:52:10 2010 UTC revision 2972 by artak, Thu Mar 4 01:33:38 2010 UTC
# Line 14  Line 14 
14  \section{Steady-state Heat Refraction}  \section{Steady-state Heat Refraction}
15  \label{STEADY-STATE HEAT REFRACTION}  \label{STEADY-STATE HEAT REFRACTION}
16    
17  In this chapter we will learn how to handle more complex geometries.  In this chapter we show how to handle more complex geometries.
18    
19  Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex models. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modeling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in \esc to define properties like material constants and source locations.  Steady-state heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in \esc to define properties like material constants and source locations.
20    
21  We will introduce in two stages. First we will look at a very simple geometry. In fact,  we will solve  We proceed in this chapter by first looking at a very simple geometry. In fact, we solve
22  the steady-state heat equation over a rectangular domain. This is not a very interesting problem but in a second  the steady-state heat equation over a rectangular domain. This is not a very interesting problem but then we extend our script by introducing an internal, curved interface.
 step we will extend our script by introducing an internal, curved interface.  
23    
24  \begin{figure}[ht]  \begin{figure}[ht]
25  \centerline{\includegraphics[width=4.in]{figures/pycadrec}}  \centerline{\includegraphics[width=4.in]{figures/pycadrec}}
# Line 28  step we will extend our script by introd Line 27  step we will extend our script by introd
27  \label{fig:pycad rec}  \label{fig:pycad rec}
28  \end{figure}  \end{figure}
29    
30  \section{Example 4: Creating the model with \pycad}  \section{Example 4: Creating the domain with \pycad}
31  \sslist{example04a.py}  \sslist{example04a.py}
32    
33  We will modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: firstly we will look the steady state  We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look the steady state
34  case with slightly modified boundary conditions and secondly we will use a more flexible tool  case with slightly modified boundary conditions and use a more flexible tool
35  to generate to generate the geometry. Lets look at the geometry first.  to generate to generate the geometry. Lets look at the geometry first.
36    
37  We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and  We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and
38  insulation to both sides we want to calculate the steady-state temperature distribution.  insulation to both sides we want to calculate the steady-state temperature distribution.
39    
40  In \pycad there are a few primary constructors that build upon each other to define domains and boundaries;  In \pycad there are a few primary constructors that build upon each other to define domains and boundaries;
41  the ones we will use are:  the ones we use are:
42  \begin{python}  \begin{python}
43  from esys.pycad import *  from esys.pycad import *
44  Point() #Create a point in space.  Point() #Create a point in space.
# Line 88  an instance of the \verb|Design| class w Line 87  an instance of the \verb|Design| class w
87  from esys.pycad.gmsh import Design  from esys.pycad.gmsh import Design
88  d=Design(dim=2, element_size=200*m)  d=Design(dim=2, element_size=200*m)
89  \end{python}  \end{python}
90  The argument \verb|dim| defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verb|element_size| defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be choose with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m  The argument \verb|dim| defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verb|element_size| defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be chosen with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m
91  and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So will end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily.  and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So we end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily.
92  We can easily add our domain \verb|rec| to the \verb|Design|;  We can easily add our domain \verb|rec| to the \verb|Design|;
93  \begin{python}  \begin{python}
94  d.addItem(rec)  d.addItem(rec)
95  \end{python}  \end{python}
96  We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this  We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this
97  but \pycad offers a more convenient technique called tagging. With this techniques items in the domain are  but \pycad offers a more convenient technique called tagging. With this technique items in the domain are
98  named using the \verb|PropertySet| class. We can then later use this name to set values secifically for  named using the \verb|PropertySet| class. We can then later use this name to set values secifically for
99  those sample points located on the named items. Here we name the bottom face of the  those sample points located on the named items. Here we name the bottom face of the
100  domain where we will set the heat influx:  domain where we will set the heat influx:
# Line 117  way~\footnote{An alternative are the \te Line 116  way~\footnote{An alternative are the \te
116  # write domain to an text file  # write domain to an text file
117  domain.write("example04.fly")  domain.write("example04.fly")
118  \end{python}  \end{python}
119  and then for recovery in another script;  and then for reading in another script;
120  \begin{python}  \begin{python}
121  # read domain from text file  # read domain from text file
122  from esys.finley import ReadMesh  from esys.finley import ReadMesh
# Line 132  domain =ReadMesh("example04.fly") Line 131  domain =ReadMesh("example04.fly")
131    
132  Before we discuss how we solve the PDE for the  Before we discuss how we solve the PDE for the
133  problem it is useful to present two additional options of the \verb|Design|.  problem it is useful to present two additional options of the \verb|Design|.
134  They allow the user accessing the script which is used to by \gmsh to generate the mesh as well as  They allow the user accessing the script which is used by \gmsh to generate the mesh as well as
135  the mesh as it has been generated by \gmsh. This is done by setting specific names for these files:  the mesh as it has been generated by \gmsh. This is done by setting specific names for these files:
136  \begin{python}  \begin{python}
137  d.setScriptFileName("example04.geo")  d.setScriptFileName("example04.geo")
# Line 140  d.setMeshFileName("example04.msh") Line 139  d.setMeshFileName("example04.msh")
139  \end{python}  \end{python}
140  Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and  Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and
141  the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage.  the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage.
142  Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesher can be visualise with \gmsh from the command line using  Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using
143  \begin{verbatim}  \begin{verbatim}
144  gmsh example04.geo  # show geometry  gmsh example04.geo  # show geometry
145  gmsh example04.msh  # show mesh  gmsh example04.msh  # show mesh
# Line 156  The mesh is shown in \reffig{fig:pycad r Line 155  The mesh is shown in \reffig{fig:pycad r
155  \section{The Steady-state Heat Equation}  \section{The Steady-state Heat Equation}
156  \sslist{example04b.py, cblib}  \sslist{example04b.py, cblib}
157    
158  Steady state is reached in the temperature is not changing in time. So to calculate the  Steady state is reached in the temperature when it is not changing in time. So to calculate the
159  the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero;  the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero;
160  \begin{equation}\label{eqn:Tform nabla steady}  \begin{equation}\label{eqn:Tform nabla steady}
161  -\nabla \cdot \kappa \nabla T = q\hackscore H  -\nabla \cdot \kappa \nabla T = q\hackscore H
162  \end{equation}  \end{equation}
163  This PDE is easier then the PDE in \refEq{eqn:hdgenf2} we needed to solve in each time step. We can just drop  This PDE is easier to solve, than the PDE in \refEq{eqn:hdgenf2} in each time step. We can just drop
164  the \verb|D| term;  the \verb|D| term;
165  \begin{python}  \begin{python}
166  mypde=LinearPDE(domain)  mypde=LinearPDE(domain)
# Line 173  already discussed how this constraint is Line 172  already discussed how this constraint is
172  x=Solution(domain).getX()  x=Solution(domain).getX()
173  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
174  \end{python}  \end{python}
175  Notice that we have use the \verb|sup| function to calculate the maximum $y$ coordinator of the relevant sample points.  Notice that we use the \verb|sup| function to calculate the maximum of $y$ coordinates of the relevant sample points.
176    
177  In all cases so far we have assumed that the domain is insulated which translates  In all cases so far we have assumed that the domain is insulated which translates
178  into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling  into a zero normal flux $-n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling
# Line 186  maintain insulation at the left and righ Line 185  maintain insulation at the left and righ
185  where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero  where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero
186  for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face.  for locations on the left or right face of the domain while it has the value \verb|qin| at the bottom face.
187  Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here.  Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here.
188  One could use the masking techniques we have already extensively used to define $q\hackscore{S}$ in \esc  We could define $q\hackscore{S}$ by using the masking techniques demonstrated earlier. The tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise
 but the tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise  
189  constant functions such as $q\hackscore{S}$. You need to remember now that we  constant functions such as $q\hackscore{S}$. You need to remember now that we
190  marked the bottom face with the name \verb|linebottom| when we defined the domain.  marked the bottom face with the name \verb|linebottom| when we defined the domain.
191  We can use this now to create $q\hackscore{S}$;  We can use this now to create $q\hackscore{S}$;
# Line 235  via the \verb|toRegGrid| function in the Line 233  via the \verb|toRegGrid| function in the
233  Our heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependent upon the material properties and shape.  Our heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependent upon the material properties and shape.
234    
235    
236  We will modify the example of the previous section by subdividing the block into two parts. The curve  We modify the example of the previous section by subdividing the block into two parts. The curve
237  separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points  separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points
238  used to define the curve may be measurement but for simplicity we assume here that there coordinates are  used to define the curve may be measurement but for simplicity we assume here that there coordinates are
239  known in analytic form.  known in analytic form.
# Line 244  There are two modes available in our exa Line 242  There are two modes available in our exa
242    
243  It is now possible to start defining our domain and boundaries.  It is now possible to start defining our domain and boundaries.
244    
245  The curve defining our clinal structure is located in approximately the middle of the domain and has a sinusoidal shape. We define the curve by generating points at discrete intervals; $51$ in this case, and then create a smooth curve through the points using the \verb Spline()  function.  The curve defining our clinal structure is located approximately in the middle of the domain and has a sinusoidal shape. We define the curve by generating points at discrete intervals; $51$ in this case, and then create a smooth curve through the points using the \verb Spline()  function.
246  \begin{python}  \begin{python}
247  # Material Boundary  # Material Boundary
248  x=[ Point(i*dsp\  x=[ Point(i*dsp\
# Line 277  This reverse spline option unfortunately Line 275  This reverse spline option unfortunately
275  #clockwise check  #clockwise check
276  bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))  bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1))
277  \end{python}  \end{python}
278  The last few steps in creating the model take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc.  The last few steps in creating the model are: take the previously defined domain and sub-domain points and generate a mesh that can be imported into \esc.
279  To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim  and a specified number of finite elements that need to be applied to the domain \verb element_size  . It then becomes a simple task of adding the sub-domains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the sub-domain properties in the solution script. This is done using the \verb PropertySet()  function. The geometry and mesh are then saved so the \esc domain can be created.  To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim  and a specified number of finite elements that need to be applied to the domain \verb element_size  . It then becomes a simple task of adding the sub-domains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the sub-domain properties in the solution script. This is done using the \verb PropertySet()  function. The geometry and mesh are then saved so the \esc domain can be created.
280  \begin{python}  \begin{python}
281  # Create a Design which can make the mesh  # Create a Design which can make the mesh
# Line 315  No further changes are required to the P Line 313  No further changes are required to the P
313  \sslist{example05b.py and cblib.py}  \sslist{example05b.py and cblib.py}
314    
315  We want to investigate the profile of the data of the last example.  We want to investigate the profile of the data of the last example.
316  We are particularly interested in the deepth profile of the heat flux which is  We are particularly interested in the depth profile of the heat flux which is
317  the second component of $-\kappa \nabla T$. We will extend the script developed in the  the second component of $-\kappa \nabla T$. We extend the script developed in the
318  previous section to show how for instance vertical profile can be plotted.  previous section to show how for instance vertical profile can be plotted.
319    
320  The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and  The first important information is that \esc assumes that $-\kappa \nabla T$ is not smooth and
# Line 324  will in fact represent the values at num Line 322  will in fact represent the values at num
322  the flux is the product of the piecewise constant function $\kappa$ and  the flux is the product of the piecewise constant function $\kappa$ and
323  the gradient of the temperature $T$ which has a kink across the rock interface.  the gradient of the temperature $T$ which has a kink across the rock interface.
324  Before plotting this function we need to smooth it using the  Before plotting this function we need to smooth it using the
325  \verb|Projector()| factory;  \verb|Projector()| class;
326  \begin{python}  \begin{python}
327  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector
328  proj=Projector(domain)  proj=Projector(domain)
# Line 337  to a regular $200\times 200$ grid; Line 335  to a regular $200\times 200$ grid;
335  \begin{python}  \begin{python}
336  xiq,yiq,ziq = toRegGrid(qu[1],200,200)  xiq,yiq,ziq = toRegGrid(qu[1],200,200)
337  \end{python}  \end{python}
338  using the \verb|toRegGrid| function from the cookbook library which we are already using for the contour plot.  using the \verb|toRegGrid| function from the cookbook library which we are using for the contour plot.
339  At return \verb|ziq[j,i]| is the value of vertical heat flux at point  At return \verb|ziq[j,i]| is the value of vertical heat flux at point
340  (\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by  (\verb|xiq[i]|,\verb|yiq[j]|). We can easily create deep profiles now by
341  plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script  plotting slices \verb|ziq[:,i]| over \verb|yiq|. The following script
# Line 360  This process can be repeated for various Line 358  This process can be repeated for various
358    
359  \section{Arrow plots in \mpl}  \section{Arrow plots in \mpl}
360  \sslist{example05c.py and cblib.py}  \sslist{example05c.py and cblib.py}
361  We would like to visualize the distribution of the flux $-\kappa \nabla T$ over the domain  We would like to visualise the distribution of the flux $-\kappa \nabla T$ over the domain
362  and produce a plot like shown in \reffig{fig:hr001qumodel}.  and produce a plot like shown in \reffig{fig:hr001qumodel}.
363  The plot puts together three components. A contour plot of the temperature,  The plot puts together three components. A contour plot of the temperature,
364  a colored representation of the two sub-domains where colour represents the thermal conductivity  a colored representation of the two sub-domains where colour represents the thermal conductivity
Legend:
Removed from v.2971  
changed lines
  Added in v.2972
  ViewVC Help
Powered by ViewVC 1.1.26