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

Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: application/x-tex
File size: 12338 byte(s)
Make everyone sad by touching all the files


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