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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3427 - (hide annotations)
Thu Jan 6 00:11:45 2011 UTC (8 years, 9 months ago) by ahallam
File MIME type: text/x-python
File size: 3331 byte(s)
Missing script.
1 ahallam 3427
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.finley import Rectangle # This imports the rectangle domain function from finley
36     import os, sys #This package is necessary to handle saving our data.
37     from math import pi, sqrt, sin, cos
38    
39     from esys.escript.pdetools import Projector
40     from cblib import toRegGrid
41    
42     import matplotlib
43     matplotlib.use('agg') #It's just here for automated testing
44    
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     #Domain related.
56     mx = 1000*m #meters - model length
57     my = -250*m #meters - model depth
58     ndx = 200 # mesh steps in x direction
59     ndy = 200 # mesh steps in y direction - one dimension means one element
60     #PDE related
61     res1=500.0
62     res2=10.0
63     res3=10000.0
64     #con=1/res
65     cur=1000.
66    
67     ################################################ESTABLISHING PARAMETERS
68     #the folder to put our outputs in, leave blank "" for script path
69     save_path= os.path.join("data","example11")
70     #ensure the dir exists
71     mkDir(save_path)
72    
73     ####################################################DOMAIN CONSTRUCTION
74     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
75     x=Solution(domain).getX()
76    
77     kro=kronecker(domain)
78     source1=[3.*mx/8.,0]; source2=[5.*mx/8.,0]
79    
80     c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)
81     c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)
82     sourceg=cur*(c1-c2)
83    
84     res=res1*wherePositive(x[1]-my/3)+res2*whereNegative(x[1]-my/3)*wherePositive(x[1]-my*2/3)+res3*whereNegative(x[1]-my*2/3)
85     con=1/res
86     q=whereZero(x[1]-my)+whereZero(x[0])+whereZero(x[0]-mx)
87     ###############################################ESCRIPT PDE CONSTRUCTION
88    
89     mypde=LinearPDE(domain)
90     mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)
91     #mypde.setSymmetryOn()
92     sol=mypde.getSolution()
93    
94     # Save the output to file.
95     saveVTK(os.path.join(save_path,"ex11b.vtu"),\
96     source=sourceg,\
97     res_pot=sol,\
98     res=res,\
99     curden=-con*grad(sol),\
100     efield=-grad(sol))

  ViewVC Help
Powered by ViewVC 1.1.26