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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months ago) by jfenwick
File MIME type: text/x-python
File size: 3929 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 jfenwick 3981 ##############################################################################
2 ahallam 3135 #
3 jfenwick 6651 # Copyright (c) 2009-2018 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ahallam 3135 #
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 3135 #
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 3135
17 jfenwick 6651 __copyright__="""Copyright (c) 2009-2018 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3135 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 3135 __url__="https://launchpad.net/escript-finley"
23    
24     """
25     Author: Antony Hallam antony.hallam@uqconnect.edu.au
26     """
27    
28     ############################################################FILE HEADER
29     # example10a.py
30 ahallam 3195 # Model of gravitational Potential for a Gravity Well.
31 ahallam 3135
32     #######################################################EXTERNAL MODULES
33     # To solve the problem it is necessary to import the modules we require.
34     from esys.escript import * # This imports everything from the escript library
35     from esys.escript.unitsSI import *
36     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
37 caltinay 3346 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
38 ahallam 3135 import os, sys #This package is necessary to handle saving our data.
39    
40 ahallam 3389 from esys.escript.pdetools import Projector, Locator
41     import numpy as np
42    
43 sshaw 5288 try:
44     # This imports the rectangle domain function
45     from esys.finley import Brick
46     HAVE_FINLEY = True
47     except ImportError:
48     print("Finley module not available")
49     HAVE_FINLEY = False
50 ahallam 3135 ########################################################MPI WORLD CHECK
51     if getMPISizeWorld() > 1:
52     import sys
53 jfenwick 3892 print("This example will not run in an MPI world.")
54 ahallam 3135 sys.exit(0)
55    
56 sshaw 5288 if HAVE_FINLEY:
57     #################################################ESTABLISHING VARIABLES
58     #Domain related.
59     mx = 500*m #meters - model length
60     my = 500*m #meters - model width
61     mz = -4000*m
62     ndx = 15 # mesh steps in x direction
63     ndy = 15 # mesh steps in y direction - one dimension means one element
64     ndz = 10
65     #PDE related
66     rho=100.0
67     rholoc=[0,0,-0]
68     G=6.67300*10E-11
69 ahallam 3135
70 sshaw 5288 ################################################ESTABLISHING PARAMETERS
71     #the folder to put our outputs in, leave blank "" for script path
72     save_path= os.path.join("data","example10")
73     #ensure the dir exists
74     mkDir(save_path)
75 ahallam 3135
76 sshaw 5288 ####################################################DOMAIN CONSTRUCTION
77     domain = Brick(l0=mx,l1=my,n0=ndx, n1=ndy,l2=mz,n2=ndz)
78     x=Solution(domain).getX()
79     x=x-[mx/2,my/2,mz/2]
80     domain.setX(interpolate(x, ContinuousFunction(domain)))
81     mask=wherePositive(100-length(x-rholoc))
82     rho=rho*mask
83     kro=kronecker(domain)
84 ahallam 3135
85 sshaw 5288 mass=rho*vol(domain)
86     ipot=FunctionOnBoundary(domain)
87     xb=ipot.getX()
88 ahallam 3389
89 sshaw 5288 q=whereZero(x[2]-inf(x[2]))
90     ###############################################ESCRIPT PDE CONSTRUCTION
91 ahallam 3135
92 sshaw 5288 mypde=LinearPDE(domain)
93     mypde.setValue(A=kro,Y=4.*3.1415*G*rho,q=q,r=0)
94     mypde.setSymmetryOn()
95     sol=mypde.getSolution()
96     saveVTK(os.path.join(save_path,"ex10b.vtu"),\
97     grav_pot=sol,\
98     g_field=-grad(sol),\
99     g_fieldz=-grad(sol)*[0,0,1],\
100     gz=length(-grad(sol)*[0,0,1]))
101 ahallam 3135
102 sshaw 5288 ################################################MODEL SIZE SAMPLING
103     sampler=[]
104     for i in range(-250,250,1):
105     sampler.append([i,0,250])
106 ahallam 3389
107 sshaw 5288 sample=[] # array to hold values
108     rec=Locator(domain,sampler) #location to record
109     psol=rec.getValue(sol)
110     np.savetxt(os.path.join(save_path,"example10b_%04d.asc"%mx),psol)
111 ahallam 3389

  ViewVC Help
Powered by ViewVC 1.1.26