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

Diff of /trunk/doc/examples/cookbook/example05b.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1 
1  from __future__ import division  from __future__ import division, print_function
 from __future__ import print_function  
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2014 by University of Queensland  # Copyright (c) 2009-2014 by University of Queensland
# Line 36  import matplotlib Line 35  import matplotlib
35  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
36  from esys.pycad import * #domain constructor  from esys.pycad import * #domain constructor
37  from esys.pycad.gmsh import Design #Finite Element meshing package  from esys.pycad.gmsh import Design #Finite Element meshing package
 from esys.finley import MakeDomain #Converter for escript  
38  import os #file path tool  import os #file path tool
39  from math import * # math package  from math import * # math package
40  from esys.escript import *  from esys.escript import *
# Line 46  from esys.escript.pdetools import Projec Line 44  from esys.escript.pdetools import Projec
44  from cblib import toRegGrid  from cblib import toRegGrid
45  import pylab as pl #Plotting package  import pylab as pl #Plotting package
46    
47    try:
48        # This imports the rectangle domain function
49        from esys.finley import MakeDomain#Converter for escript
50        HAVE_FINLEY = True
51    except ImportError:
52        print("Finley module not available")
53        HAVE_FINLEY = False
54  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
55  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
56          import sys          import sys
57          print("This example will not run in an MPI world.")          print("This example will not run in an MPI world.")
58          sys.exit(0)          sys.exit(0)
59    if HAVE_FINLEY:
60  #################################################ESTABLISHING VARIABLES      #################################################ESTABLISHING VARIABLES
61  #set modal to 1 for a syncline or -1 for an anticline structural      #set modal to 1 for a syncline or -1 for an anticline structural
62  #configuration      #configuration
63  modal=-1      modal=-1
64    
65  # the folder to put our outputs in, leave blank "" for script path -      # the folder to put our outputs in, leave blank "" for script path -
66  # note this folder path must exist to work      # note this folder path must exist to work
67  save_path= os.path.join("data","example05")      save_path= os.path.join("data","example05")
68  mkDir(save_path)      mkDir(save_path)
69    
70  ################################################ESTABLISHING PARAMETERS      ################################################ESTABLISHING PARAMETERS
71  #Model Parameters      #Model Parameters
72  width=5000.0*m   #width of model      width=5000.0*m   #width of model
73  depth=-6000.0*m  #depth of model      depth=-6000.0*m  #depth of model
74  Ttop=20*K       # top temperature      Ttop=20*K       # top temperature
75  qin=70*Milli*W/(m*m) # bottom heat influx      qin=70*Milli*W/(m*m) # bottom heat influx
76    
77  sspl=51 #number of discrete points in spline      sspl=51 #number of discrete points in spline
78  dsp=width/(sspl-1) #dx of spline steps for width      dsp=width/(sspl-1) #dx of spline steps for width
79  dep_sp=2500.0*m #avg depth of spline      dep_sp=2500.0*m #avg depth of spline
80  h_sp=1500.0*m #heigh of spline      h_sp=1500.0*m #heigh of spline
81  orit=-1.0 #orientation of spline 1.0=>up -1.0=>down      orit=-1.0 #orientation of spline 1.0=>up -1.0=>down
82    
83  ####################################################DOMAIN CONSTRUCTION      ####################################################DOMAIN CONSTRUCTION
84  # Domain Corners      # Domain Corners
85  p0=Point(0.0,      0.0, 0.0)      p0=Point(0.0,      0.0, 0.0)
86  p1=Point(0.0,    depth, 0.0)      p1=Point(0.0,    depth, 0.0)
87  p2=Point(width, depth, 0.0)      p2=Point(width, depth, 0.0)
88  p3=Point(width,   0.0, 0.0)      p3=Point(width,   0.0, 0.0)
89    
90  # Generate Material Boundary      # Generate Material Boundary
91  x=[ Point(i*dsp\      x=[ Point(i*dsp\
92      ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\          ,-dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\
93       for i in range(0,sspl)\           for i in range(0,sspl)\
94    ]        ]
95  mysp = Spline(*tuple(x))      mysp = Spline(*tuple(x))
96  # Start and end of material boundary.      # Start and end of material boundary.
97  x1=mysp.getStartPoint()      x1=mysp.getStartPoint()
98  x2=mysp.getEndPoint()      x2=mysp.getEndPoint()
99    
100  #  Create TOP BLOCK      #  Create TOP BLOCK
101  # lines      # lines
102  tbl1=Line(p0,x1)      tbl1=Line(p0,x1)
103  tbl2=mysp      tbl2=mysp
104  tbl3=Line(x2,p3)      tbl3=Line(x2,p3)
105  l30=Line(p3, p0)      l30=Line(p3, p0)
106  # curve      # curve
107  tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)      tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30)
108  # surface      # surface
109  tblock = PlaneSurface(tblockloop)      tblock = PlaneSurface(tblockloop)
110    
111    
112  # Create BOTTOM BLOCK      # Create BOTTOM BLOCK
113  # lines      # lines
114  bbl1=Line(x1,p1)      bbl1=Line(x1,p1)
115  bbl3=Line(p2,x2)      bbl3=Line(p2,x2)
116  bbl4=-mysp      bbl4=-mysp
117  l12=Line(p1, p2)      l12=Line(p1, p2)
118  # curve      # curve
119  bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)      bblockloop = CurveLoop(bbl1,l12,bbl3,bbl4)
120  # surface      # surface
121  bblock = PlaneSurface(bblockloop)      bblock = PlaneSurface(bblockloop)
122    
123  ################################################CREATE MESH FOR ESCRIPT      ################################################CREATE MESH FOR ESCRIPT
124  # Create a Design which can make the mesh      # Create a Design which can make the mesh
125  d=Design(dim=2, element_size=200)      d=Design(dim=2, element_size=200)
126  # Add the subdomains and flux boundaries.      # Add the subdomains and flux boundaries.
127  d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\      d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
128                                       PropertySet("linebottom",l12))                                           PropertySet("linebottom",l12))
129  # Create the geometry, mesh and Escript domain      # Create the geometry, mesh and Escript domain
130  d.setScriptFileName(os.path.join(save_path,"example05.geo"))      d.setScriptFileName(os.path.join(save_path,"example05.geo"))
131  d.setMeshFileName(os.path.join(save_path,"example05.msh"))      d.setMeshFileName(os.path.join(save_path,"example05.msh"))
132  domain=MakeDomain(d, optimizeLabeling=True)      domain=MakeDomain(d, optimizeLabeling=True)
133  print("Domain has been generated ...")      print("Domain has been generated ...")
134  ##############################################################SOLVE PDE      ##############################################################SOLVE PDE
135  mypde=LinearPDE(domain)      mypde=LinearPDE(domain)
136  mypde.getSolverOptions().setVerbosityOn()      mypde.getSolverOptions().setVerbosityOn()
137  mypde.setSymmetryOn()      mypde.setSymmetryOn()
138  kappa=Scalar(0,Function(domain))      kappa=Scalar(0,Function(domain))
139  kappa.setTaggedValue("top",2.0*W/m/K)      kappa.setTaggedValue("top",2.0*W/m/K)
140  kappa.setTaggedValue("bottom",4.0*W/m/K)      kappa.setTaggedValue("bottom",4.0*W/m/K)
141  mypde.setValue(A=kappa*kronecker(domain))      mypde.setValue(A=kappa*kronecker(domain))
142  x=Solution(domain).getX()      x=Solution(domain).getX()
143  mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)      mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
144  qS=Scalar(0,FunctionOnBoundary(domain))      qS=Scalar(0,FunctionOnBoundary(domain))
145  qS.setTaggedValue("linebottom",qin)      qS.setTaggedValue("linebottom",qin)
146  mypde.setValue(y=qS)      mypde.setValue(y=qS)
147  print("PDE has been generated ...")      print("PDE has been generated ...")
148  ###########################################################GET SOLUTION      ###########################################################GET SOLUTION
149  T=mypde.getSolution()      T=mypde.getSolution()
150  print("PDE has been solved  ...")      print("PDE has been solved  ...")
151    
152  ##################################################REGRIDDING & PLOTTING      ##################################################REGRIDDING & PLOTTING
153  xi, yi, zi = toRegGrid(T, nx=50, ny=50)      xi, yi, zi = toRegGrid(T, nx=50, ny=50)
154  pl.matplotlib.pyplot.autumn()      pl.matplotlib.pyplot.autumn()
155  pl.contourf(xi,yi,zi,10)      pl.contourf(xi,yi,zi,10)
156  pl.xlabel("Horizontal Displacement (m)")      pl.xlabel("Horizontal Displacement (m)")
157  pl.ylabel("Depth (m)")      pl.ylabel("Depth (m)")
158  pl.savefig(os.path.join(save_path,"Tcontour.png"))      pl.savefig(os.path.join(save_path,"Tcontour.png"))
159  print("Solution has been plotted  ...")      print("Solution has been plotted  ...")
160  ##########################################################VISUALISATION      ##########################################################VISUALISATION
161  # calculate gradient of solution for quiver plot      # calculate gradient of solution for quiver plot
162  #Projector is used to smooth the data.      #Projector is used to smooth the data.
163  proj=Projector(domain)      proj=Projector(domain)
164  #move data to a regular grid for plotting      #move data to a regular grid for plotting
165  xi,yi,zi = toRegGrid(T,200,200)      xi,yi,zi = toRegGrid(T,200,200)
166  cut=int(len(xi)//2)      cut=int(len(xi)//2)
167  pl.clf()      pl.clf()
168  pl.plot(zi[:,cut],yi)      pl.plot(zi[:,cut],yi)
169  pl.title("Temperature Depth Profile")      pl.title("Temperature Depth Profile")
170  pl.xlabel("Temperature (K)")      pl.xlabel("Temperature (K)")
171  pl.ylabel("Depth (m)")      pl.ylabel("Depth (m)")
172  pl.savefig(os.path.join(save_path,"tdp.png"))      pl.savefig(os.path.join(save_path,"tdp.png"))
173  pl.clf()      pl.clf()
174                
175  # Heat flow depth profile.      # Heat flow depth profile.
176  # grid the data.      # grid the data.
177  qu=proj(-kappa*grad(T))      qu=proj(-kappa*grad(T))
178  xiq,yiq,ziq = toRegGrid(qu[1],50,50)      xiq,yiq,ziq = toRegGrid(qu[1],50,50)
179  cut=int(len(xiq)//2)      cut=int(len(xiq)//2)
180  pl.plot(ziq[:,cut]*1000.,yiq)      pl.plot(ziq[:,cut]*1000.,yiq)
181  pl.title("Vertical Heat Flow Depth Profile")      pl.title("Vertical Heat Flow Depth Profile")
182  pl.xlabel("Heat Flow (mW/m^2)")      pl.xlabel("Heat Flow (mW/m^2)")
183  pl.ylabel("Depth (m)")      pl.ylabel("Depth (m)")
184  pl.savefig(os.path.join(save_path,"hf.png"))      pl.savefig(os.path.join(save_path,"hf.png"))
185  pl.clf()      pl.clf()
186    
187  # Temperature Gradient Depth Profile at x[50]      # Temperature Gradient Depth Profile at x[50]
188  zT=proj(-grad(T))      zT=proj(-grad(T))
189  xt,yt,zt=toRegGrid(zT[1],200,200)      xt,yt,zt=toRegGrid(zT[1],200,200)
190  cut=int(len(xt)//2)      cut=int(len(xt)//2)
191  pl.plot(zt[:,cut]*1000.,yt)      pl.plot(zt[:,cut]*1000.,yt)
192  pl.title("Vertical Temperature Gradient \n Depth Profile")      pl.title("Vertical Temperature Gradient \n Depth Profile")
193  pl.xlabel("Temperature gradient (K/Km)")      pl.xlabel("Temperature gradient (K/Km)")
194  pl.ylabel("Depth (m)")      pl.ylabel("Depth (m)")
195  pl.savefig(os.path.join(save_path,"tgdp.png"))      pl.savefig(os.path.join(save_path,"tgdp.png"))
196  pl.clf()      pl.clf()
197    
198  # Thermal Conditions Depth Profile          # Thermal Conditions Depth Profile    
199  xk,yk,zk = toRegGrid(proj(kappa),200,200)      xk,yk,zk = toRegGrid(proj(kappa),200,200)
200  cut=int(len(xk)//2)      cut=int(len(xk)//2)
201  pl.plot(zk[:,cut],yk)      pl.plot(zk[:,cut],yk)
202  pl.title("Thermal Conductivity Depth Profile")      pl.title("Thermal Conductivity Depth Profile")
203  pl.xlabel("Conductivity (W/K/m)")      pl.xlabel("Conductivity (W/K/m)")
204  pl.ylabel("Depth (m)")      pl.ylabel("Depth (m)")
205  pl.axis([1,5,-6000,0])      pl.axis([1,5,-6000,0])
206  pl.savefig(os.path.join(save_path,"tcdp.png"))      pl.savefig(os.path.join(save_path,"tcdp.png"))
207  pl.clf()      pl.clf()
208  print("vertical profiles created ...")      print("vertical profiles created ...")

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26