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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26