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

Annotation of /trunk/doc/examples/cookbook/example04b.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2948 - (hide annotations)
Thu Feb 25 04:54:30 2010 UTC (9 years, 1 month ago) by gross
Original Path: trunk/doc/examples/cookbook/simple_solver001.py
File MIME type: text/x-python
File size: 3695 byte(s)
a new almost completed version of the cookbook
1 gross 2948
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     # simple_mesher001.py
28     # Create either a 2D mesh for a rectangle model using pycad meshing
29     # tools.
30     #
31     #######################################################EXTERNAL MODULES
32     from esys.pycad import * #domain constructor
33     from esys.pycad.gmsh import Design #Finite Element meshing package
34     from esys.finley import MakeDomain #Converter for escript
35     from esys.escript import *
36     from esys.escript.unitsSI import *
37     from esys.escript.linearPDEs import LinearPDE
38     import matplotlib
39     matplotlib.use('agg') #It's just here for automated testing
40     import pylab as pl #Plotting package
41     from cblib import toRegGrid
42     import os
43    
44     ########################################################MPI WORLD CHECK
45     if getMPISizeWorld() > 1:
46     import sys
47     print "This example will not run in an MPI world."
48     sys.exit(0)
49    
50     # make sure path exists
51     save_path= os.path.join("data","simpleheat")
52     mkDir(save_path)
53    
54     ################################################ESTABLISHING PARAMETERS
55     #Model Parameters
56     width=5000.0*m #width of model
57     depth=-6000.0*m #depth of model
58     kappa=2.0*W/m/K # watts/m.Kthermal conductivity
59     Ttop=20*K # top temperature
60     qin=70*Milli*W/(m*m) # bottom heat influx
61    
62     ####################################################DOMAIN CONSTRUCTION
63     # Domain Corners
64     p0=Point(0.0, 0.0, 0.0)
65     p1=Point(0.0, depth, 0.0)
66     p2=Point(width, depth, 0.0)
67     p3=Point(width, 0.0, 0.0)
68     # Join corners in anti-clockwise manner.
69     l01=Line(p0, p1)
70     l12=Line(p1, p2)
71     l23=Line(p2, p3)
72     l30=Line(p3, p0)
73     # Join line segments to create domain boundary.
74     c=CurveLoop(l01, l12, l23, l30)
75     # surface
76     rec = PlaneSurface(c)
77    
78     #############################################EXPORTING MESH FOR ESCRIPT
79     # Create a Design which can make the mesh
80     d=Design(dim=2, element_size=200*m)
81     # Add the subdomains and flux boundaries.
82     d.addItems(rec, PropertySet("linebottom",l12))
83     #############################################MAKE THE FINLEY DOMAIN
84     domain=MakeDomain(d, optimizeLabeling=True)
85     print "Domain has been generated ..."
86     ############################################# solve PDE
87     mypde=LinearPDE(domain)
88     mypde.getSolverOptions().setVerbosityOn()
89     mypde.setSymmetryOn()
90     mypde.setValue(A=kappa*kronecker(domain))
91     x=Solution(domain).getX()
92     mypde.setValue(q=whereZero(x[1]-sup(x[1])),r=Ttop)
93     qS=Scalar(0,FunctionOnBoundary(domain))
94     qS.setTaggedValue("linebottom",qin)
95     mypde.setValue(y=-qS)
96     print "PDE has been generated ..."
97     ###########################################################GET SOLUTION
98     T=mypde.getSolution()
99     print "PDE has been solved ..."
100     ###########################################################
101     xi, yi, zi = toRegGrid(T, nx=50, ny=50)
102     pl.matplotlib.pyplot.autumn()
103     pl.contourf(xi,yi,zi,10)
104     pl.xlabel("Horizontal Displacement (m)")
105     pl.ylabel("Depth (m)")
106     pl.savefig(os.path.join(save_path,"simpleheat.png"))
107     print "Solution has been plotted ..."

  ViewVC Help
Powered by ViewVC 1.1.26