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

Annotation of /trunk/doc/examples/cookbook/example10d.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3389 - (hide annotations)
Wed Dec 1 23:03:35 2010 UTC (8 years, 10 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 ahallam 3389
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