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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3099 - (hide annotations)
Sun Aug 22 04:35:15 2010 UTC (8 years, 6 months ago) by ahallam
File MIME type: text/x-python
File size: 3247 byte(s)
Gravitational Potential Example- Needs some tuning.
1 ahallam 3099
2     ########################################################
3     #
4     # Copyright (c) 2009 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 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     # example10a.py
28     # Model of gravitational Potential.
29    
30     #######################################################EXTERNAL MODULES
31     # To solve the problem it is necessary to import the modules we require.
32     from esys.escript import * # This imports everything from the escript library
33     from esys.escript.unitsSI import *
34     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
35     from esys.escript.linearPDEs import Poisson # This defines LinearPDE as LinearPDE
36     from esys.finley import Rectangle # This imports the rectangle domain function from finley
37     #For interactive use, you can comment out the next two lines
38     import matplotlib
39     matplotlib.use('agg') #It's just here for automated testing
40     import pylab as pl #Plotting package.
41     import numpy as np #Array package.
42     import os, sys #This package is necessary to handle saving our data.
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     #################################################ESTABLISHING VARIABLES
51     #Domain related.
52     mx = 500*m #meters - model length
53     my = 500*m #meters - model width
54     ndx = 500 # mesh steps in x direction
55     ndy = 500 # mesh steps in y direction - one dimension means one element
56     #PDE related
57     rho=100.0
58     rholoc=[250,250]
59     G=6.67300*10E-11
60    
61     ################################################ESTABLISHING PARAMETERS
62     #the folder to put our outputs in, leave blank "" for script path
63     save_path= os.path.join("data","example10a")
64     #ensure the dir exists
65     mkDir(save_path)
66    
67     ####################################################DOMAIN CONSTRUCTION
68     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
69     x=Solution(domain).getX()
70     #mask=whereZero(x[0]-rholoc[0])*whereZero(x[1]-rholoc[1])
71     mask=wherePositive(100-length(x-rholoc))
72     #rho=Scalar(rho,ContinuousFunction(domain))
73     rho=rho*mask
74     y=Scalar(0.0,FunctionOnBoundary(domain))
75     kro=kronecker(domain)
76    
77     q=whereZero(x[0])+\
78     whereZero(x[1])+\
79     whereZero(x[0]-sup(x[0]))+\
80     whereZero(x[1]-sup(x[1]))
81     ###############################################ESCRIPT PDE CONSTRUCTION
82    
83     mypde=LinearPDE(domain)
84     mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0.)
85     #mypde.setSymmetryOn()
86     #mypde=Poisson(domain)
87     #mypde.setValue(f=rho*4.*3.1415*G,q=q)
88     sol=mypde.getSolution()
89     saveVTK("ex10a.vtu",field_strength=sol)#,field=-grad(sol))

  ViewVC Help
Powered by ViewVC 1.1.26