/[escript]/trunk/doc/examples/cookbook/example06.py
ViewVC logotype

Contents of /trunk/doc/examples/cookbook/example06.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4853 - (show annotations)
Wed Apr 9 05:41:57 2014 UTC (5 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 5593 byte(s)
Added import division and changed a few / to //

1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2009-2014 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 __copyright__="""Copyright (c) 2009-2014 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 """
26 Author: Antony Hallam antony.hallam@uqconnect.edu.au
27 """
28 ############################################################FILE HEADER
29 # example06.py
30 # Create a 3 block fault and overburden style model in 2d using pycad
31 # meshing tools.
32
33 #######################################################EXTERNAL MODULES
34 import matplotlib
35 matplotlib.use('agg') #It's just here for automated testing
36 from esys.pycad import *
37 from esys.pycad.gmsh import Design
38 from esys.finley import MakeDomain
39 from esys.escript import *
40 import numpy as np
41 import pylab as pl #Plotting package
42 from cblib import toRegGrid, subsample
43 from esys.escript.unitsSI import *
44 from esys.escript.linearPDEs import LinearPDE
45 import os, sys
46
47 ########################################################MPI WORLD CHECK
48 if getMPISizeWorld() > 1:
49 import sys
50 print("This example will not run in an MPI world.")
51 sys.exit(0)
52 #################################################ESTABLISHING VARIABLES
53 # where to put output files
54 save_path= os.path.join("data","example06")
55 mkDir(save_path)
56
57 Ttop=20*Celsius # temperature at the top
58 qin=300.*Milli*W/(m*m) #our heat source temperature is now zero
59
60 ################################################ESTABLISHING PARAMETERS
61 # Overall Domain
62 p0=Point(0.0, 0.0, 0.0)
63 p1=Point(0.0, -6000.0, 0.0)
64 p2=Point(5000.0, -6000.0, 0.0)
65 p3=Point(5000.0, 0.0, 0.0)
66
67 ####################################################DOMAIN CONSTRUCTION
68 l01=Line(p0, p1)
69 l12=Line(p1, p2)
70 l23=Line(p2, p3)
71 l30=Line(p3, p0)
72
73 # Generate Material Boundary
74 p4=Point(0.0, -2400.0, 0.0)
75 p5=Point(2000.0, -2400.0, 0.0)
76 p6=Point(3000.0, -6000.0, 0.0)
77 p7=Point(5000.0, -2400.0, 0.0)
78
79 # Create TOP BLOCK
80 tbl1=Line(p0,p4)
81 tbl2=Line(p4,p5)
82 tbl3=Line(p5,p7)
83 tbl4=Line(p7,p3)
84 tblockloop = CurveLoop(tbl1,tbl2,tbl3,tbl4,l30)
85 tblock = PlaneSurface(tblockloop)
86
87 # Create BOTTOM BLOCK LEFT
88 bbll1=Line(p4,p1)
89 bbll2=Line(p1,p6)
90 bbll3=Line(p6,p5)
91 bbll4=-tbl2
92 bblockloopl = CurveLoop(bbll1,bbll2,bbll3,bbll4)
93 bblockl = PlaneSurface(bblockloopl)
94
95 # Create BOTTOM BLOCK RIGHT
96 bbrl1=Line(p6,p2)
97 bbrl2=Line(p2,p7)
98 bbrl3=-tbl3
99 bbrl4=-bbll3
100 bblockloopr = CurveLoop(bbrl1,bbrl2,bbrl3,bbrl4)
101 bblockr = PlaneSurface(bblockloopr)
102
103 #############################################EXPORTING MESH FOR ESCRIPT
104 # Create a Design which can make the mesh
105 d=Design(dim=2, element_size=200)
106 # Add the subdomains and flux boundaries.
107 d.addItems(PropertySet("top",tblock),\
108 PropertySet("bottomleft",bblockl),\
109 PropertySet("bottomright",bblockr),\
110 PropertySet("linebottom",bbll2, bbrl1))
111
112 # Create the geometry, mesh and Escript domain
113 d.setScriptFileName(os.path.join(save_path,"example06.geo"))
114
115 d.setMeshFileName(os.path.join(save_path,"example06.msh"))
116 domain=MakeDomain(d)
117 print("Domain has been generated ...")
118
119 # set up kappa (thermal conductivity across domain) using tags
120 kappa=Scalar(0,Function(domain))
121 kappa.setTaggedValue("top",2.0)
122 kappa.setTaggedValue("bottomleft",10.0)
123 kappa.setTaggedValue("bottomright",6.0)
124 ##############################################################SOLVE PDE
125 mypde=LinearPDE(domain)
126 mypde.getSolverOptions().setVerbosityOn()
127 mypde.setSymmetryOn()
128 mypde.setValue(A=kappa*kronecker(domain))
129 x=Solution(domain).getX()
130 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
131 qS=Scalar(0,FunctionOnBoundary(domain))
132 qS.setTaggedValue("linebottom",qin)
133 mypde.setValue(y=qS)
134 print("PDE has been generated ...")
135 ###########################################################GET SOLUTION
136 T=mypde.getSolution()
137 print("PDE has been solved ...")
138 ###############################################################PLOTTING
139 # show temperature:
140 xi, yi, zi = toRegGrid(T, nx=50, ny=50)
141 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
142 pl.clabel(CS, inline=1, fontsize=8)
143 # show sub domains:
144 tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
145 pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
146 bpgr=np.array([p.getCoordinates() for p in bblockloopr.getPolygon() ])
147 pl.fill(bpgr[:,0],bpgr[:,1],'orange',label='6 W/m/k',zorder=-1000)
148 bpgl=np.array([p.getCoordinates() for p in bblockloopl.getPolygon() ])
149 pl.fill(bpgl[:,0],bpgl[:,1],'red',label='10 W/m/k',zorder=-1000)
150 # show flux:
151 xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
152 pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
153 # create plot
154 pl.title("Heat Refraction across an anisotropic structure\n with flux vectors.")
155 pl.xlabel("Horizontal Displacement (m)")
156 pl.ylabel("Depth (m)")
157 pl.legend()
158 pl.savefig(os.path.join(save_path,"flux.png"))
159 print("Flux has been plotted ...")

  ViewVC Help
Powered by ViewVC 1.1.26