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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (hide annotations)
Tue Apr 10 08:57:23 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 3551 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



1 ahallam 3426
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 jfenwick 3892 print("This example will not run in an MPI world.")
53 ahallam 3426 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