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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 4957 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 # example10d.py
28 # Model of gravitational Potential for a gravity POLE.
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 from esys.weipa import saveVTK # This imports the VTK file saver from weipa
37 import os, sys #This package is necessary to handle saving our data.
38 from math import pi, sqrt, sin, cos
39
40 from esys.escript.pdetools import Projector, Locator
41 from esys.finley import ReadGmsh
42
43 import matplotlib
44 matplotlib.use('agg') #Its just here for automated testing
45
46 import pylab as pl #Plotting package
47 import numpy as np
48 from mpl_toolkits.mplot3d import Axes3D
49
50 ########################################################MPI WORLD CHECK
51 if getMPISizeWorld() > 1:
52 import sys
53 print("This example will not run in an MPI world.")
54 sys.exit(0)
55
56 #################################################ESTABLISHING VARIABLES
57 #Domain related.
58 mx = 10000*m #meters - model length
59 my = 10000*m #meters - model width
60 #PDE related
61 rho=10.0
62 rholoc=[mx/2.,my/2.]
63 G=6.67300*10E-11
64
65 R=10
66 z=50.
67
68 ################################################ESTABLISHING PARAMETERS
69 #the folder to put our outputs in, leave blank "" for script path
70 save_path= os.path.join("data","example10")
71 #ensure the dir exists
72 mkDir(save_path)
73
74 #####################################################ANALYTIC SOLUTION
75 def analytic_gz(x,z,R,drho):
76 G=6.67300*10E-11
77 return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
78
79 sol_angz=[]
80 sol_anx=[]
81 for x in range(int(-mx/20),int(mx/20),10):
82 sol_angz.append(analytic_gz(x,z,R,rho))
83 sol_anx.append(x+mx/2)
84
85 ##############INVERSION
86 def gzpot(p, y, x, *args):
87 #rho, rhox, rhoy, R = p
88 rhox=args[0]/2.; rhoy=args[1]/2.
89 rho, R, z =p
90 #Domain related.
91 mx = args[0]; my = args[1];
92 #PDE related
93 G=6.67300*10E-11
94
95 #DOMAIN CONSTRUCTION
96 domain=ReadGmsh('data/example10m/example10m.msh',2)
97 domx=Solution(domain).getX()
98 mask=wherePositive(R-length(domx-rholoc))
99 rhoe=rho*mask
100 kro=kronecker(domain)
101
102 q=whereZero(domx[1]-my)+whereZero(domx[1])+whereZero(domx[0])+whereZero(domx[0]-mx)
103 #ESCRIPT PDE CONSTRUCTION
104 mypde=LinearPDE(domain)
105 mypde.setValue(A=kro,Y=4.*np.pi*G*rhoe,q=q,r=0.0)
106 mypde.setSymmetryOn()
107 sol=mypde.getSolution()
108
109 g_field=grad(sol) #The graviational accelleration g.
110 g_fieldz=g_field*[0,1] #The vertical component of the g field.
111 gz=length(g_fieldz) #The magnitude of the vertical component.
112
113 #MODEL SIZE SAMPLING
114 sol_escgz=[]
115 sol_escx=[]
116 for i in range(0,len(x)):
117 sol_escgz.append([x[i],rhoy+z])
118
119 sample=[] # array to hold values
120 rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
121 psol=rec.getValue(gz)
122
123 err = np.sum((np.array(y) - np.array(psol))**2.)
124 print("Lsup= ",Lsup(np.array(psol)-np.array(sol_angz))/Lsup(np.array(psol)))
125 return err
126
127 #Initial Guess
128 #guess=[400,mx/4,my/4,50]
129 guess=[15.,20.]
130
131 from scipy.optimize import leastsq
132 #plsq = leastsq(gzpot, guess, args=(sol_angz, sol_anx, mx, my, ndx, ndy),maxfev=20)
133 #print plsq
134
135 objf=[]
136 x=np.arange(1,20)
137 y=np.arange(1,20)
138 z=np.arange(40,60)
139 fig=pl.figure(figsize=(5,5))
140
141 for p in x:
142 objf.append(gzpot([p,10.,50.],sol_angz,sol_anx, mx, my))
143 sp=fig.add_subplot(311)
144 sp.plot(x,objf)
145 sp.set_title("Variable RHO")
146 objf=[]
147 for R in y:
148 objf.append(gzpot([10.,R,50.],sol_angz,sol_anx, mx, my))
149 sp=fig.add_subplot(312)
150 sp.plot(y,objf)
151 sp.set_title("Variable Radius")
152
153 objf=[]
154 for D in z:
155 objf.append(gzpot([10.,10.,D],sol_angz,sol_anx, mx, my))
156 sp=fig.add_subplot(313)
157 sp.plot(z,objf)
158 sp.set_title("Variable Depth")
159
160 fig.savefig("ex10e_objf.pdf",dpi=150)
161
162 #ob=np.array(objf)
163 #X,Y=pl.meshgrid(x,y)
164 #fig=pl.figure()
165 #ax=Axes3D(fig)
166 #ax.plot_surface(X,Y,ob)
167
168 #pl.show()
169

  ViewVC Help
Powered by ViewVC 1.1.26