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

  ViewVC Help
Powered by ViewVC 1.1.26