/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 6508 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26