/[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 5288 - (show annotations)
Tue Dec 2 23:18:40 2014 UTC (4 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 5669 byte(s)
fixing tests for cases where required domains not built
1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2009-2014 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 __copyright__="""Copyright (c) 2009-2014 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 """
26 Author: Antony Hallam antony.hallam@uqconnect.edu.au
27 """
28
29 ############################################################FILE HEADER
30 # example10d.py
31 # Model of gravitational Potential for a gravity POLE.
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') #Its 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 esys.finley import ReadGmsh
47
48 import pylab as pl #Plotting package
49 import numpy as np
50 from mpl_toolkits.mplot3d import Axes3D
51
52 try:
53 # This imports the rectangle domain function
54 from esys.finley import MakeDomain
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 #Domain related.
68 mx = 10000*m #meters - model length
69 my = 10000*m #meters - model width
70 #PDE related
71 rho=10.0
72 rholoc=[mx/2.,my/2.]
73 G=6.67300*10E-11
74
75 R=10
76 z=50.
77
78 ################################################ESTABLISHING PARAMETERS
79 #the folder to put our outputs in, leave blank "" for script path
80 save_path= os.path.join("data","example10")
81 #ensure the dir exists
82 mkDir(save_path)
83
84 #####################################################ANALYTIC SOLUTION
85 def analytic_gz(x,z,R,drho):
86 G=6.67300*10E-11
87 return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
88
89 sol_angz=[]
90 sol_anx=[]
91 for x in range(int(-mx/20),int(mx/20),10):
92 sol_angz.append(analytic_gz(x,z,R,rho))
93 sol_anx.append(x+mx/2)
94
95 ##############INVERSION
96 def gzpot(p, y, x, *args):
97 #rho, rhox, rhoy, R = p
98 rhox=args[0]/2.; rhoy=args[1]/2.
99 rho, R, z =p
100 #Domain related.
101 mx = args[0]; my = args[1];
102 #PDE related
103 G=6.67300*10E-11
104
105 #DOMAIN CONSTRUCTION
106 domain=ReadGmsh('data/example10m/example10m.msh',2)
107 domx=Solution(domain).getX()
108 mask=wherePositive(R-length(domx-rholoc))
109 rhoe=rho*mask
110 kro=kronecker(domain)
111
112 q=whereZero(domx[1]-my)+whereZero(domx[1])+whereZero(domx[0])+whereZero(domx[0]-mx)
113 #ESCRIPT PDE CONSTRUCTION
114 mypde=LinearPDE(domain)
115 mypde.setValue(A=kro,Y=4.*np.pi*G*rhoe,q=q,r=0.0)
116 mypde.setSymmetryOn()
117 sol=mypde.getSolution()
118
119 g_field=grad(sol) #The graviational accelleration g.
120 g_fieldz=g_field*[0,1] #The vertical component of the g field.
121 gz=length(g_fieldz) #The magnitude of the vertical component.
122
123 #MODEL SIZE SAMPLING
124 sol_escgz=[]
125 sol_escx=[]
126 for i in range(0,len(x)):
127 sol_escgz.append([x[i],rhoy+z])
128
129 sample=[] # array to hold values
130 rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
131 psol=rec.getValue(gz)
132
133 err = np.sum((np.array(y) - np.array(psol))**2.)
134 print("Lsup= ",Lsup(np.array(psol)-np.array(sol_angz))/Lsup(np.array(psol)))
135 return err
136
137 #Initial Guess
138 #guess=[400,mx/4,my/4,50]
139 guess=[15.,20.]
140
141 from scipy.optimize import leastsq
142 #plsq = leastsq(gzpot, guess, args=(sol_angz, sol_anx, mx, my, ndx, ndy),maxfev=20)
143 #print plsq
144
145 objf=[]
146 x=np.arange(1,20)
147 y=np.arange(1,20)
148 z=np.arange(40,60)
149 fig=pl.figure(figsize=(5,5))
150
151 for p in x:
152 objf.append(gzpot([p,10.,50.],sol_angz,sol_anx, mx, my))
153 sp=fig.add_subplot(311)
154 sp.plot(x,objf)
155 sp.set_title("Variable RHO")
156 objf=[]
157 for R in y:
158 objf.append(gzpot([10.,R,50.],sol_angz,sol_anx, mx, my))
159 sp=fig.add_subplot(312)
160 sp.plot(y,objf)
161 sp.set_title("Variable Radius")
162
163 objf=[]
164 for D in z:
165 objf.append(gzpot([10.,10.,D],sol_angz,sol_anx, mx, my))
166 sp=fig.add_subplot(313)
167 sp.plot(z,objf)
168 sp.set_title("Variable Depth")
169
170 fig.savefig("ex10e_objf.pdf",dpi=150)
171
172 #ob=np.array(objf)
173 #X,Y=pl.meshgrid(x,y)
174 #fig=pl.figure()
175 #ax=Axes3D(fig)
176 #ax.plot_surface(X,Y,ob)
177
178 #pl.show()
179

  ViewVC Help
Powered by ViewVC 1.1.26