/[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 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4861 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 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 for an infinite strike horizontal
29 # cylinder. This example tests the boundary condition criterion.
30
31 #######################################################EXTERNAL MODULES
32 # To solve the problem it is necessary to import the modules we require.
33 from esys.escript import * # This imports everything from the escript library
34 from esys.escript.unitsSI import *
35 from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
36 from esys.finley import Rectangle # This imports the rectangle domain function from finley
37 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
38 import os, sys #This package is necessary to handle saving our data.
39 from math import pi, sqrt, sin, cos
40
41 from esys.escript.pdetools import Projector, Locator
42 from cblib import toRegGrid
43 from esys.finley import ReadGmsh
44
45 import matplotlib
46 matplotlib.use('agg') #It's just here for automated testing
47
48 import pylab as pl #Plotting package
49 import numpy as np
50
51 ########################################################MPI WORLD CHECK
52 if getMPISizeWorld() > 1:
53 import sys
54 print("This example will not run in an MPI world.")
55 sys.exit(0)
56
57 #################################################ESTABLISHING VARIABLES
58 #PDE related
59 rho=200.0 #the density contrast of the source
60 G=6.67300*10E-11 #gravitational constant
61 R=100. #radius of the source body
62 ################################################ESTABLISHING PARAMETERS
63 #the folder to put our outputs in, leave blank "" for script path
64 save_path= os.path.join("data","example10")
65 mesh_path = "data/example10m"
66 #ensure the dir exists
67 mkDir(save_path)
68
69 ####################################################DOMAIN CONSTRUCTION
70 #Geometric and material property related variables.
71 domain=ReadGmsh(os.path.join(mesh_path,'example10m.msh'),2) # create the domain
72 x=domain.getX()
73 #Domain related dimensions from the imported file.
74 mx = Lsup(x)*m #meters - model length
75 my = Lsup(x)*m #meters - model width
76 rholoc=[mx/2,my/2]
77
78 mask=wherePositive(R-length(x-rholoc))
79 rhoe=rho*mask
80 kro=kronecker(domain)
81
82 #Define the analytic solution.
83 def analytic_gz(x,z,R,drho):
84 G=6.67300*10E-11
85 return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
86
87 #Define the boundary conditions.
88 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[0]-mx)
89 ###############################################ESCRIPT PDE CONSTRUCTION
90 mypde=LinearPDE(domain)
91 mypde.setValue(A=kro,Y=4.*3.1415*G*rhoe,q=q,r=0.)
92 mypde.setSymmetryOn()
93 sol=mypde.getSolution()
94
95 g_field=grad(sol) #The gravitational acceleration g.
96 g_fieldz=g_field*[0,1] #The vertical component of the g field.
97 gz=length(g_fieldz) #The magnitude of the vertical component.
98 # Save the output to file.
99 saveVTK(os.path.join(save_path,"ex10d.vtu"),\
100 grav_pot=sol,g_field=g_field,g_fieldz=g_fieldz,gz=gz,rho=rhoe)
101
102 ################################################MODEL SIZE SAMPLING
103 smoother = Projector(domain) #Function smoother.
104 z=1000. #Distance of profile from source.
105 sol_angz=[] #Array for analytic gz
106 sol_anx=[] #Array for x
107 #calculate analytic gz and x location.
108 for x in range(int(-mx/2.),int(mx/2.),10):
109 sol_angz.append(analytic_gz(x,z,R,rho))
110 sol_anx.append(x+mx/2)
111 #save analytic solution
112 pl.savetxt(os.path.join(save_path,"ex10d_as.asc"),sol_angz)
113
114 sol_escgz=[] #Array for escript solution for gz
115 #calculate the location of the profile in the domain
116 for i in range(0,len(sol_anx)):
117 sol_escgz.append([sol_anx[i],my/2.-z])
118
119 rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
120 psol=rec.getValue(smoother(gz)) #extract the solution
121 #save X and escript solution profile to file
122 pl.savetxt(os.path.join(save_path,"ex10d_%05d.asc"%mx),psol)
123 pl.savetxt(os.path.join(save_path,"ex10d_x%05d.asc"%mx),sol_anx)
124
125 #plot the pofile and analytic solution using example10q.py
126

  ViewVC Help
Powered by ViewVC 1.1.26