/[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 3389 - (show annotations)
Wed Dec 1 23:03:35 2010 UTC (8 years, 8 months ago) by ahallam
File MIME type: text/x-python
File size: 4100 byte(s)
Some changes to the cookbook examples. Starting Inversion experimentation.
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 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 cblib import toRegGrid
42
43 import matplotlib
44 matplotlib.use('agg') #It's just here for automated testing
45
46 import pylab as pl #Plotting package
47 import numpy as np
48
49 ########################################################MPI WORLD CHECK
50 if getMPISizeWorld() > 1:
51 import sys
52 print "This example will not run in an MPI world."
53 sys.exit(0)
54
55 #################################################ESTABLISHING VARIABLES
56 #Domain related.
57 mx = 4000*m #meters - model length
58 my = -4000*m #meters - model width
59 ndx = 400 # mesh steps in x direction
60 ndy = 400 # mesh steps in y direction - one dimension means one element
61 #PDE related
62 rho=200.0
63 rholoc=[mx/2.,my/2.]
64 G=6.67300*10E-11
65
66 R=10
67 z=50.
68
69 ################################################ESTABLISHING PARAMETERS
70 #the folder to put our outputs in, leave blank "" for script path
71 save_path= os.path.join("data","example10")
72 #ensure the dir exists
73 mkDir(save_path)
74
75 #####################################################ANALYTIC SOLUTION
76 def analytic_gz(x,z,R,drho):
77 G=6.67300*10E-11
78 return G*2*np.pi*R*R*drho*(z/(x*x+z*z))
79
80 sol_angz=[]
81 sol_anx=[]
82 for x in range(-mx/2,mx/2,10):
83 sol_angz.append(analytic_gz(x,z,R,rho))
84 sol_anx.append(x+mx/2)
85
86 ##############INVERSION
87 def gzpot(p, y, x, *args):
88 rho, rhox, rhoy, R = p
89
90 #Domain related.
91 mx = args[0]; my = args[1]; ndx = args[2]; ndy = args[3]
92 #PDE related
93 G=6.67300*10E-11
94
95 #DOMAIN CONSTRUCTION
96 domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
97 x=Solution(domain).getX()
98 mask=wherePositive(R-length(x-rholoc))
99 rhoe=rho*mask
100 kro=kronecker(domain)
101
102 q=whereZero(x[1]-my)+whereZero(x[1])+whereZero(x[0])+whereZero(x[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],my/2.+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 = y - psol
124
125 #Initial Guess
126 guess=[400,mx/4,my/4,50]
127
128 from scipy.optimize import leastsq
129 plsq = leastsq(gzpot, guess, args=(sol_angz, sol_anx, mx, my, ndx, ndy))
130

  ViewVC Help
Powered by ViewVC 1.1.26