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

  ViewVC Help
Powered by ViewVC 1.1.26