/[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 6892 - (hide annotations)
Tue Oct 15 04:55:57 2019 UTC (8 days, 19 hours ago) by aellery
File MIME type: text/x-python
File size: 3868 byte(s)
- Updated the jupyter dockerfile
- Fixed the problem that John was having with the examples


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

  ViewVC Help
Powered by ViewVC 1.1.26