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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 months ago) by sshaw
File MIME type: text/x-python
File size: 5422 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2009-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __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 # Model of gravitational Potential for an infinite strike horizontal
31 # cylinder. This example tests the boundary condition criterion.
32
33 #######################################################EXTERNAL MODULES
34 # To solve the problem it is necessary to import the modules we require.
35 import matplotlib
36 matplotlib.use('agg') #It's just here for automated testing
37
38 from esys.escript import * # This imports everything from the escript library
39 from esys.escript.unitsSI import *
40 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
41 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
42 import os, sys #This package is necessary to handle saving our data.
43 from math import pi, sqrt, sin, cos
44
45 from esys.escript.pdetools import Projector, Locator
46 from cblib import toRegGrid
47 from esys.finley import ReadGmsh
48
49 import pylab as pl #Plotting package
50 import numpy as np
51
52 try:
53 # This imports the rectangle domain function
54 from esys.finley import Rectangle
55 HAVE_FINLEY = True
56 except ImportError:
57 print("Finley module not available")
58 HAVE_FINLEY = False
59 ########################################################MPI WORLD CHECK
60 if getMPISizeWorld() > 1:
61 import sys
62 print("This example will not run in an MPI world.")
63 sys.exit(0)
64
65 if HAVE_FINLEY:
66 #################################################ESTABLISHING VARIABLES
67 #PDE related
68 rho=200.0 #the density contrast of the source
69 G=6.67300*10E-11 #gravitational constant
70 R=100. #radius of the source body
71 ################################################ESTABLISHING PARAMETERS
72 #the folder to put our outputs in, leave blank "" for script path
73 save_path= os.path.join("data","example10")
74 mesh_path = "data/example10m"
75 #ensure the dir exists
76 mkDir(save_path)
77
78 ####################################################DOMAIN CONSTRUCTION
79 #Geometric and material property related variables.
80 domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
81 x=domain.getX()
82 #Domain related dimensions from the imported file.
83 mx = Lsup(x)*m #meters - model length
84 my = Lsup(x)*m #meters - model width
85 rholoc=[mx/2,my/2]
86
87 mask=wherePositive(R-length(x-rholoc))
88 rhoe=rho*mask
89 kro=kronecker(domain)
90
91 #Define the analytic solution.
92 def analytic_gz(x,z,R,drho):
93 G=6.67300*10E-11
94 return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
95
96 #Define the boundary conditions.
97 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
98 ###############################################ESCRIPT PDE CONSTRUCTION
99 mypde=LinearPDE(domain)
100 mypde.setValue(A=kro,Y=4.*3.1415*G*rhoe,q=q,r=0.)
101 mypde.setSymmetryOn()
102 sol=mypde.getSolution()
103
104 g_field=grad(sol) #The gravitational acceleration g.
105 g_fieldz=g_field*[0,1] #The vertical component of the g field.
106 gz=length(g_fieldz) #The magnitude of the vertical component.
107 # Save the output to file.
108 saveVTK(os.path.join(save_path,"ex10d.vtu"),\
109 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz,rho=rhoe)
110
111 ################################################MODEL SIZE SAMPLING
112 smoother = Projector(domain) #Function smoother.
113 z=1000. #Distance of profile from source.
114 sol_angz=[] #Array for analytic gz
115 sol_anx=[] #Array for x
116 #calculate analytic gz and x location.
117 for x in range(int(-mx/2.),int(mx/2.),10):
118 sol_angz.append(analytic_gz(x,z,R,rho))
119 sol_anx.append(x+mx/2)
120 #save analytic solution
121 pl.savetxt(os.path.join(save_path,"ex10d_as.asc"),sol_angz)
122
123 sol_escgz=[] #Array for escript solution for gz
124 #calculate the location of the profile in the domain
125 for i in range(0,len(sol_anx)):
126 sol_escgz.append([sol_anx[i],my/2.-z])
127
128 rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
129 psol=rec.getValue(smoother(gz)) #extract the solution
130 #save X and escript solution profile to file
131 pl.savetxt(os.path.join(save_path,"ex10d_%05d.asc"%mx),psol)
132 pl.savetxt(os.path.join(save_path,"ex10d_x%05d.asc"%mx),sol_anx)
133
134 #plot the pofile and analytic solution using example10q.py
135

  ViewVC Help
Powered by ViewVC 1.1.26