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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 5606 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1
2 ########################################################
3 #
4 # Copyright (c) 2009-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 __copyright__="""Copyright (c) 2009-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 """
23 Author: Antony Hallam antony.hallam@uqconnect.edu.au
24 """
25
26 ############################################################FILE HEADER
27 # example05c.py
28 # Create either a 2D syncline or anticline model using pycad meshing
29 # tools. The model steady state heat solution. How to create arrow or
30 # quiver plots.
31
32 #######################################################EXTERNAL MODULES
33 import matplotlib
34 matplotlib.use('agg') #It's just here for automated testing
35 from esys.pycad import * #domain constructor
36 from esys.pycad.gmsh import Design #Finite Element meshing package
37 from esys.finley import MakeDomain #Converter for escript
38 import os #file path tool
39 from math import * # math package
40 from esys.escript import *
41 from esys.escript.unitsSI import *
42 from esys.escript.linearPDEs import LinearPDE
43 from esys.escript.pdetools import Projector
44 from cblib import toRegGrid, subsample
45 import pylab as pl #Plotting package
46 import numpy as np
47
48 ########################################################MPI WORLD CHECK
49 if getMPISizeWorld() > 1:
50 import sys
51 print("This example will not run in an MPI world.")
52 sys.exit(0)
53
54 #################################################ESTABLISHING VARIABLES
55 #set modal to 1 for a syncline or -1 for an anticline structural
56 #configuration
57 modal=-1
58
59 # the folder to put our outputs in, leave blank "" for script path -
60 # note this folder path must exist to work
61 save_path= os.path.join("data","example05")
62 mkDir(save_path)
63
64 ################################################ESTABLISHING PARAMETERS
65 #Model Parameters
66 width=5000.0*m #width of model
67 depth=-6000.0*m #depth of model
68 Ttop=20*K # top temperature
69 qin=70*Milli*W/(m*m) # bottom heat influx
70
71 sspl=51 #number of discrete points in spline
72 dsp=width/(sspl-1) #dx of spline steps for width
73 dep_sp=2500.0*m #avg depth of spline
74 h_sp=1500.0*m #heigh of spline
75 orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
76
77 ####################################################DOMAIN CONSTRUCTION
78 # Domain Corners
79 p0=Point(0.0, 0.0, 0.0)
80 p1=Point(0.0, depth, 0.0)
81 p2=Point(width, depth, 0.0)
82 p3=Point(width, 0.0, 0.0)
83
84 # Generate Material Boundary
85 x=[ Point(i*dsp\
86 ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
87 for i in range(0,sspl)\
88 ]
89 mysp = Spline(*tuple(x))
90 # Start and end of material boundary.
91 x1=mysp.getStartPoint()
92 x2=mysp.getEndPoint()
93
94 # Create TOP BLOCK
95 # lines
96 tbl1=Line(p0,x1)
97 tbl2=mysp
98 tbl3=Line(x2,p3)
99 l30=Line(p3, p0)
100 # curve
101 tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
102 # surface
103 tblock = PlaneSurface(tblockloop)
104 # Create BOTTOM BLOCK
105 # lines
106 bbl1=Line(x1,p1)
107 bbl3=Line(p2,x2)
108 bbl4=-mysp
109 l12=Line(p1, p2)
110 # curve
111 bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
112
113 # surface
114 bblock = PlaneSurface(bblockloop)
115
116 ################################################CREATE MESH FOR ESCRIPT
117 # Create a Design which can make the mesh
118 d=Design(dim=2, element_size=200)
119 # Add the subdomains and flux boundaries.
120 d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
121 PropertySet("linebottom",l12))
122 # Create the geometry, mesh and Escript domain
123 d.setScriptFileName(os.path.join(save_path,"example05.geo"))
124 d.setMeshFileName(os.path.join(save_path,"example05.msh"))
125 domain=MakeDomain(d, optimizeLabeling=True)
126 print("Domain has been generated ...")
127 ##############################################################SOLVE PDE
128 mypde=LinearPDE(domain)
129 mypde.getSolverOptions().setVerbosityOn()
130 mypde.setSymmetryOn()
131 kappa=Scalar(0,Function(domain))
132 kappa.setTaggedValue("top",2.0*W/m/K)
133 kappa.setTaggedValue("bottom",4.0*W/m/K)
134 mypde.setValue(A=kappa*kronecker(domain))
135 x=Solution(domain).getX()
136 mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
137 qS=Scalar(0,FunctionOnBoundary(domain))
138 qS.setTaggedValue("linebottom",qin)
139 mypde.setValue(y=qS)
140 print("PDE has been generated ...")
141 ###########################################################GET SOLUTION
142 T=mypde.getSolution()
143 print("PDE has been solved ...")
144 ###############################################################PLOTTING
145 # show temperature:
146 xi, yi, zi = toRegGrid(T, nx=50, ny=50)
147 CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k')
148 pl.clabel(CS, inline=1, fontsize=8)
149 # show sub domains:
150 tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ])
151 pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=-1000)
152 bpg=np.array([p.getCoordinates() for p in bblockloop.getPolygon() ])
153 pl.fill(bpg[:,0],bpg[:,1],'red',label='4 W/m/k',zorder=-1000)
154 # show flux:
155 xflux, flux=subsample(-kappa*grad(T), nx=20, ny=20)
156 pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white")
157 # create plot
158 pl.title("Heat Refraction across a clinal structure\n with heat flux.")
159 pl.xlabel("Horizontal Displacement (m)")
160 pl.ylabel("Depth (m)")
161 pl.legend()
162 pl.savefig(os.path.join(save_path,"flux.png"))
163 print("Flux has been plotted ...")

  ViewVC Help
Powered by ViewVC 1.1.26