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

Contents of /trunk/doc/cookbook/example04.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3153 - (show annotations)
Sun Sep 5 01:31:49 2010 UTC (10 years, 4 months ago) by ahallam
File MIME type: application/x-tex
File size: 11923 byte(s)
Updates to cookbook. Fixes to SWP and new Gravitational Potential
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % Earth Systems Science Computational Center (ESSCC)
6 % http://www.uq.edu.au/esscc
7 %
8 % Primary Business: Queensland, Australia
9 % Licensed under the Open Software License version 3.0
10 % http://www.opensource.org/licenses/osl-3.0.php
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14 \section{Steady-state Heat Refraction}
15 \label{STEADY-STATE HEAT REFRACTION}
16
17 In this chapter we demonstrate how to handle more complex geometries.
18
19 Steady-state heat refraction will give us an opportunity to investigate some of
20 the richer features that the \esc package has to offer. One of these is \pycad .
21 The advantage of using \pycad is that it offers an easy method for developing
22 and manipulating complex domains. In conjunction with \gmsh we can generate
23 finite element meshes that conform to our domain's shape providing accurate
24 modelling of interfaces and boundaries. Another useful function of \pycad is
25 that we can tag specific areas of our domain with labels as we construct them.
26 These labels can then be used in \esc to define properties like material
27 constants and source locations.
28
29 We proceed in this chapter by first looking at a very simple geometry. Whilst a
30 simple rectangular domain is not very interesting the example is elaborated upon
31 later by introducing an internal curved interface.
32
33 \section{Example 4: Creating the Domain with \pycad}
34 \label{example4}
35 \sslist{example04a.py}
36 We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look at the
37 steady state case with slightly modified boundary conditions and use a more
38 flexible tool to generate the geometry. Let us look at the geometry first.
39
40 We want to define a rectangular domain of width $5 km$ and depth $6 km$ below
41 the surface of the Earth. The domain is subject to a few conditions. The
42 temperature is known at the surface and the basement has a known heat flux. Each
43 side of the domain is insulated and the aim is to calculate the final
44 temperature distribution.
45
46 In \pycad there are a few primary constructors that build upon each other to
47 define domains and boundaries. The ones we use are:
48 \begin{python}
49 from esys.pycad import *
50 Point() #Create a point in space.
51 Line() #Creates a line from a number of points.
52 CurveLoop() #Creates a closed loop from a number of lines.
53 PlaneSurface() #Creates a surface based on a CurveLoop
54 \end{python}
55 So to construct our domain as shown in \reffig{fig:pycad rec}, we first need to
56 create the corner points. From the corner points we build the four edges of the
57 rectangle. The four edges then form a closed loop which defines our domain as a
58 surface. We start by inputting the variables we need to construct the model.
59 \begin{python}
60 width=5000.0*m #width of model
61 depth=-6000.0*m #depth of model
62 \end{python}
63
64 \begin{figure}[ht]
65 \centerline{\includegraphics[width=4.in]{figures/pycadrec}}
66 \caption{Example 4: Rectangular Domain for \pycad}
67 \label{fig:pycad rec}
68 \end{figure}
69
70 The variables are then used to construct the four corners of our domain, which
71 from the origin has the dimensions of $5000$ meters width and $-6000$ meters
72 depth. This is done with the \verb|Point()| function which accepts x, y and z
73 coordinates. Our domain is in two dimensions so z should always be zero.
74 \begin{python}
75 # Overall Domain
76 p0=Point(0.0, 0.0, 0.0)
77 p1=Point(0.0, depth, 0.0)
78 p2=Point(width, depth, 0.0)
79 p3=Point(width, 0.0, 0.0)
80 \end{python}
81 Now lines are defined using our points. This forms a rectangle around our
82 domain:
83 \begin{python}
84 l01=Line(p0, p1)
85 l12=Line(p1, p2)
86 l23=Line(p2, p3)
87 l30=Line(p3, p0)
88 \end{python}
89 Note that lines have a direction. These lines form the basis for our domain
90 boundary, which is a closed loop.
91 \begin{python}
92 c=CurveLoop(l01, l12, l23, l30)
93 \end{python}
94 Be careful to define the curved loop in an \textbf{anti-clockwise} manner
95 otherwise the meshing algorithm may fail.
96 Finally we can define the domain as
97 \begin{python}
98 rec = PlaneSurface(c)
99 \end{python}
100 At this point the introduction of the curved loop seems to be unnecessary but
101 this concept plays an important role if holes are introduced.
102
103 Now we are ready to hand over the domain \verb|rec| to a mesher which
104 subdivides the domain into triangles (or tetrahedra in 3D). In our case we use
105 \gmsh. We create an instance of the \verb|Design| class which will handle the
106 interface to \gmsh:
107 \begin{python}
108 from esys.pycad.gmsh import Design
109 d=Design(dim=2, element_size=200*m)
110 \end{python}
111 The argument \verb|dim| defines the spatial dimension of the domain\footnote{If
112 \texttt{dim}=3 the rectangle would be interpreted as a surface in the three
113 dimensional space.}. The second argument \verb|element_size| defines the element
114 size which is the maximum length of a triangle edge in the mesh. The element
115 size needs to be chosen with care in order to avoid very dense meshes. If the
116 mesh is too dense, the computational time will be long but if the mesh is too
117 sparse, the modelled result will be poor. In our case with an element size of
118 $200$m and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$
119 triangles in each spatial direction. So we end up with about $30 \times 30 =
120 900$ triangles which is a size that can be handled easily.
121 The domain \verb|rec| can simply be added to the \verb|Design|;
122 \begin{python}
123 d.addItem(rec)
124 \end{python}
125 We have the plan to set a heat flux on the bottom of the domain. One can use
126 the masking technique to do this but \pycad offers a more convenient technique
127 called tagging. With this technique items in the domain are named using the
128 \verb|PropertySet| class. We can then later use this name to set values
129 specifically for those sample points located on the named items. Here we name
130 the bottom face of the domain where we will set the heat influx:
131 \begin{python}
132 ps=PropertySet("linebottom",l12))
133 d.addItem(ps)
134 \end{python}
135 Now we are ready to hand over the \verb|Design| to \FINLEY:
136 \begin{python}
137 from esys.finley import MakeDomain
138 domain=MakeDomain(d)
139 \end{python}
140 The \verb|domain| object can now be used in the same way like the return object
141 of the \verb|Rectangle| object we have used previously to generate a mesh. It
142 is common practice to separate the mesh generation from the PDE solution.
143 The main reason for this is that mesh generation can be computationally very
144 expensive in particular in 3D. So it is more efficient to generate the mesh
145 once and write it to a file. The mesh can then be read in every time a new
146 simulation is run. \FINLEY supports this in the following
147 way\footnote{An alternative is using the \texttt{dump} and \texttt{load}
148 functions. They work with a binary format and tend to be much smaller.}:
149 \begin{python}
150 # write domain to a text file
151 domain.write("example04.fly")
152 \end{python}
153 and then for reading in another script:
154 \begin{python}
155 # read domain from text file
156 from esys.finley import ReadMesh
157 domain =ReadMesh("example04.fly")
158 \end{python}
159
160 Before we discuss how to solve the PDE for this problem, it is useful to
161 present two additional options of the \verb|Design| class.
162 These allow the user to access the script which is used by \gmsh to generate
163 the mesh as well as the generated mesh itself. This is done by setting specific
164 names for these files:
165 \begin{python}
166 d.setScriptFileName("example04.geo")
167 d.setMeshFileName("example04.msh")
168 \end{python}
169 Conventionally the extension \texttt{geo} is used for the script file of the
170 \gmsh geometry and the extension \texttt{msh} for the mesh file. Normally these
171 files are deleted after usage.
172 Accessing these files can be helpful to debug the generation of more complex
173 geometries. The geometry and the mesh can be visualised from the command line
174 using
175 \begin{verbatim}
176 gmsh example04.geo # show geometry
177 gmsh example04.msh # show mesh
178 \end{verbatim}
179 The mesh is shown in \reffig{fig:pycad rec mesh}.
180
181 \begin{figure}[ht]
182 \centerline{\includegraphics[width=4.in]{figures/simplemesh}}
183 \caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}}
184 \label{fig:pycad rec mesh}
185 \end{figure}
186 \clearpage
187
188 \section{The Steady-state Heat Equation}
189 \sslist{example04b.py, cblib}
190 A temperature equilibrium or steady state is reached when the temperature
191 distribution in the model does not change with time. To calculate the steady
192 state solution the time derivative term in \refEq{eqn:Tform nabla} needs to be
193 set to zero;
194 \begin{equation}\label{eqn:Tform nabla steady}
195 -\nabla \cdot \kappa \nabla T = q\hackscore H
196 \end{equation}
197 This PDE is easier to solve than the PDE in \refEq{eqn:hdgenf2}, as no time
198 steps (iterations) are required. The \verb|D| term from \refEq{eqn:hdgenf2} is
199 simply dropped in this case.
200 \begin{python}
201 mypde=LinearPDE(domain)
202 mypde.setValue(A=kappa*kronecker(model), Y=qH)
203 \end{python}
204 The temperature at the top face of the domain is known as \verb|Ttop|~($=20 C$).
205 In \refSec{Sec:1DHDv0} we have already discussed how this constraint is added
206 to the PDE:
207 \begin{python}
208 x=Solution(domain).getX()
209 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
210 \end{python}
211 Notice that we use the \verb|sup| function to calculate the maximum of $y$
212 coordinates of the relevant sample points.
213
214 In all cases so far we have assumed that the domain is insulated which
215 translates into a zero normal flux $-n \cdot \kappa \nabla T$, see
216 \refEq{eq:hom flux}. In the modelling set-up of this chapter we want to set
217 the normal heat flux at the bottom to \verb|qin| while still maintaining
218 insulation at the left and right face. Mathematically we can express this as
219 \begin{equation}
220 -n \cdot \kappa \nabla T = q\hackscore{S}
221 \label{eq:inhom flux}
222 \end{equation}
223 where $q\hackscore{S}$ is a function of its location on the boundary. Its value
224 becomes zero for locations on the left or right face of the domain while it has
225 the value \verb|qin| at the bottom face.
226 Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we
227 prescribe the temperature here.
228 We could define $q\hackscore{S}$ by using the masking techniques demonstrated
229 earlier. The tagging mechanism provides an alternative and in many cases more
230 convenient way of defining piecewise constant functions such as
231 $q\hackscore{S}$. Recall now that the bottom face was denoted with the name
232 \verb|linebottom| when we defined the domain.
233 We can use this now to create $q\hackscore{S}$;
234 \begin{python}
235 qS=Scalar(0,FunctionOnBoundary(domain))
236 qS.setTaggedValue("linebottom",qin)
237 \end{python}
238 In the first line \verb|qS| is defined as a scalar value over the sample points
239 on the boundary of the domain. It is initialised to zero for all sample points.
240 In the second statement the values for those sample points which are located on
241 the line marked by \verb|linebottom| are set to \verb|qin|.
242
243 The Neumann boundary condition assumed by \esc has the form
244 \begin{equation}\label{NEUMAN 2b}
245 n\cdot A \cdot\nabla u = y
246 \end{equation}
247 In comparison to the version in \refEq{NEUMAN 2} we have used so far the right
248 hand side is now the new PDE coefficient $y$. As we have not specified $y$ in
249 our previous examples, \esc has assumed the value zero for $y$. A comparison of
250 \refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one needs to choose
251 $y=-q\hackscore{S}$;
252 \begin{python}
253 qS=Scalar(0,FunctionOnBoundary(domain))
254 qS.setTaggedValue("linebottom",qin)
255 mypde.setValue(y=-qS)
256 \end{python}
257 To plot the results we use the \modmpl library as shown in
258 \refSec{Sec:2DHD plot}. For convenience the interpolation of the temperature to
259 a rectangular grid for contour plotting is made available via the
260 \verb|toRegGrid| function in the \verb|cblib| module. Your result should look
261 similar to \reffig{fig:steady state heat}.
262
263 \begin{figure}[ht]
264 \centerline{\includegraphics[width=4.in]{figures/simpleheat}}
265 \caption{Example 4b: Result of simple steady state heat problem}
266 \label{fig:steady state heat}
267 \end{figure}
268 \clearpage
269

  ViewVC Help
Powered by ViewVC 1.1.26