/[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 5288 - (hide annotations)
Tue Dec 2 23:18:40 2014 UTC (4 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 3842 byte(s)
fixing tests for cases where required domains not built
1 jfenwick 4853 from __future__ import division
2 jfenwick 4848 from __future__ import print_function
3 jfenwick 3981 ##############################################################################
4 ahallam 3427 #
5 jfenwick 4657 # Copyright (c) 2009-2014 by University of Queensland
6 jfenwick 3981 # http://www.uq.edu.au
7 ahallam 3427 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
14     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 jfenwick 3981 #
16     ##############################################################################
17 ahallam 3427
18 jfenwick 4657 __copyright__="""Copyright (c) 2009-2014 by University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 ahallam 3427 Primary Business: Queensland, Australia"""
21     __license__="""Licensed under the Open Software License version 3.0
22     http://www.opensource.org/licenses/osl-3.0.php"""
23     __url__="https://launchpad.net/escript-finley"
24    
25     """
26     Author: Antony Hallam antony.hallam@uqconnect.edu.au
27     """
28    
29     ############################################################FILE HEADER
30 caltinay 4087 # example11b.py
31 ahallam 3427 # Model of gravitational Potential.
32    
33     #######################################################EXTERNAL MODULES
34     # To solve the problem it is necessary to import the modules we require.
35 caltinay 4087 import matplotlib
36     matplotlib.use('agg') #It's just here for automated testing
37 ahallam 3427 import os, sys #This package is necessary to handle saving our data.
38     from math import pi, sqrt, sin, cos
39 caltinay 4087 import numpy as np
40     import pylab as pl #Plotting package
41 ahallam 3427
42 caltinay 4087 from esys.escript import * # This imports everything from the escript library
43     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
44 ahallam 3427 from esys.escript.pdetools import Projector
45 caltinay 4087 from esys.escript.unitsSI import *
46     from esys.weipa import saveVTK
47 ahallam 3427
48 gross 3439 from cblib import toRegGrid
49    
50 sshaw 5288 try:
51     from esys.finley import Rectangle
52     HAVE_FINLEY = True
53     except ImportError:
54     print("Finley module not available")
55     HAVE_FINLEY = False
56 ahallam 3427 ########################################################MPI WORLD CHECK
57     if getMPISizeWorld() > 1:
58     import sys
59 jfenwick 3892 print("This example will not run in an MPI world.")
60 ahallam 3427 sys.exit(0)
61    
62 sshaw 5288 if HAVE_FINLEY:
63     #################################################ESTABLISHING VARIABLES
64     #Domain related.
65     mx = 1000*m #meters - model length
66     my = -250*m #meters - model depth
67     ndx = 200 # mesh steps in x direction
68     ndy = 200 # mesh steps in y direction - one dimension means one element
69     #PDE related
70     res1=500.0
71     res2=10.0
72     res3=10000.0
73     #con=1/res
74     cur=1000.
75 ahallam 3427
76 sshaw 5288 ################################################ESTABLISHING PARAMETERS
77     #the folder to put our outputs in, leave blank "" for script path
78     save_path= os.path.join("data","example11")
79     #ensure the dir exists
80     mkDir(save_path)
81 ahallam 3427
82 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
83     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
84     x=Solution(domain).getX()
85 ahallam 3427
86 sshaw 5288 kro=kronecker(domain)
87     source1=[3.*mx/8.,0]; source2=[5.*mx/8.,0]
88 ahallam 3427
89 sshaw 5288 c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)
90     c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)
91     sourceg=cur*(c1-c2)
92 ahallam 3427
93 sshaw 5288 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)
94     con=1/res
95     q=whereZero(x[1]-my)+whereZero(x[0])+whereZero(x[0]-mx)
96     ###############################################ESCRIPT PDE CONSTRUCTION
97 ahallam 3427
98 sshaw 5288 mypde=LinearPDE(domain)
99     mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)
100     #mypde.setSymmetryOn()
101     sol=mypde.getSolution()
102 ahallam 3427
103 sshaw 5288 # Save the output to file.
104     saveVTK(os.path.join(save_path,"ex11b.vtu"),\
105     source=sourceg,\
106     res_pot=sol,\
107     res=res,\
108     curden=-con*grad(sol),\
109     efield=-grad(sol))

  ViewVC Help
Powered by ViewVC 1.1.26