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

Contents of /trunk/doc/examples/cookbook/example11c.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3426 - (show annotations)
Wed Jan 5 02:51:35 2011 UTC (8 years, 3 months ago) by ahallam
File MIME type: text/x-python
File size: 3550 byte(s)
Added new example tests.
1
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 from esys.finley import ReadMesh
39
40 from esys.escript.pdetools import Projector
41 from cblib import toRegGrid
42
43 import matplotlib
44 matplotlib.use('agg') #It's just here for automated testing
45
46 import pylab as pl #Plotting package
47 import numpy as np
48
49 ########################################################MPI WORLD CHECK
50 if getMPISizeWorld() > 1:
51 import sys
52 print "This example will not run in an MPI world."
53 sys.exit(0)
54
55 #################################################ESTABLISHING VARIABLES
56 #Domain related.
57
58
59 mx = 500*m #meters - model length
60 my = 500*m #meters - model depth
61 mz = 250
62 ndx = 200 # mesh steps in x direction
63 ndy = 200 # mesh steps in y direction - one dimension means one element
64 #PDE related
65 res1=500.0
66 res2=10.0
67 res3=1000.0
68 res4=10000.0
69 #con=1/res
70 cur=1000.
71
72 ################################################ESTABLISHING PARAMETERS
73 #the folder to put our outputs in, leave blank "" for script path
74 save_path= os.path.join("data","example11")
75 #ensure the dir exists
76 mkDir(save_path)
77
78 ####################################################DOMAIN CONSTRUCTION
79 domain=ReadMesh(os.path.join(save_path,'example11lc.fly')) # create the domain
80
81 res=Scalar(0,Function(domain))
82 res.setTaggedValue("volume_0",res1)
83 res.setTaggedValue("volume_1",res2)
84 res.setTaggedValue("volume_2",res3)
85 res.setTaggedValue("volume_3",res4)
86 con=1/res
87 x=Solution(domain).getX()
88
89 kro=kronecker(domain)
90 source1=[3.*mx/8.,my/2,0]; source2=[5.*mx/8.,my/2,0]
91
92 c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)
93 c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)
94 sourceg=cur*(c1-c2)
95
96 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)+whereZero(x[2]-mz)
97 ###############################################ESCRIPT PDE CONSTRUCTION
98
99 mypde=LinearPDE(domain)
100 mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)
101 #mypde.setSymmetryOn()
102 sol=mypde.getSolution()
103
104 # Save the output to file.
105 saveVTK(os.path.join(save_path,"ex11c.vtu"),\
106 source=sourceg,\
107 res_pot=sol,\
108 res=res,\
109 curden=-con*grad(sol),\
110 abscd=length(-con*grad(sol)),\
111 efield=-grad(sol))

  ViewVC Help
Powered by ViewVC 1.1.26