/[escript]/branches/symbolic_from_3470/ripley/test/python/gravity_all.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/ripley/test/python/gravity_all.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3873 - (show annotations)
Fri Mar 16 01:12:26 2012 UTC (5 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 5184 byte(s)
a test comparing 2 aproaches for garvity
1 from esys.escript import *
2 from esys.escript.linearPDEs import LinearPDE, LinearSinglePDE
3 from esys.escript.nonlinearPDE import *
4 from esys.finley import Rectangle,Brick
5 from math import pi
6 import os
7 import esys.escript.unitsSI as U
8 #import numpy as np
9 #import scipy.optimize as so
10
11 #DIM=3
12 #NE=10
13
14 DIM=2
15 NE=80
16
17
18 class RockFeature(object):
19 def __init__(self, lx, ly, lz, x=0, y=0, depth=0, rho=0):
20 self.x=x
21 self.y=y
22 self.lx=lx
23 self.ly=ly
24 self.lz=lz
25 self.depth=depth
26 self.rho=rho
27 def getMask(self, x):
28 DIM=x.getDomain().getDim()
29 m=whereNonPositive(x[DIM-1]+self.depth) * whereNonNegative(x[DIM-1]+(self.depth+self.lz)) \
30 *whereNonNegative(x[0]-(self.x-self.lx/2)) * whereNonPositive(x[0]-(self.x+self.lx/2))
31 if DIM>2:
32 m*=whereNonNegative(x[1]-(self.y-self.ly/2)) * whereNonPositive(x[1]-(self.y+self.ly/2))
33 return m
34
35
36 H=20*U.km
37 L=60*U.km
38 L0=0
39 L1=0
40 H_earth=10.*U.km
41
42
43 rho_rock=2300*U.kg/U.m**3
44 rho_air=0.
45 feastures = [ RockFeature(lx=10.*U.km, ly=10.*U.km, lz=1.*U.km, x= L/2-10.*U.km, y=L/2-10.*U.km, depth=1.*U.km, rho=rho_rock*0.1),
46 RockFeature(lx=10.*U.km, ly=10.*U.km, lz=1.*U.km, x= L/2+10.*U.km, y=L/2+10.*U.km, depth=5.*U.km, rho=rho_rock*0.1) ]
47
48
49
50
51 # generate Domain:
52 NE_H=NE
53 NE_L=int((L/H)*NE+0.5)
54 if DIM==2:
55 domain = Rectangle(NE_L,NE_H,l0=L,l1=H)
56 x_cord=domain.getX()-[L0, H_earth]
57 else:
58 domain = Brick(NE_L,NE_L,NE_H,l0=L,l1=L,l2=H)
59 x_cord=domain.getX()-[L0, L1, H_earth]
60 LL=max(L,H)
61
62 m_psi_ref=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) + whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
63 for i in range(DIM-1):
64 m_psi_ref= m_psi_ref + whereZero(x_cord[i]-inf(x_cord[i])) + whereZero(x_cord[i]-sup(x_cord[i]))
65
66 # create test density:
67 rho_ref= 0
68 for f in feastures:
69 m=f.getMask(x_cord)
70 rho_ref = rho_ref * (1-m) + f.rho * m
71
72 #rho_ref=sup(x_cord[DIM-1])-x_cord[DIM-1] for testing!
73
74 # get the reference potential:
75 pde=LinearSinglePDE(domain)
76 pde.setValue(A=kronecker(domain), Y=4*pi*rho_ref, q=m_psi_ref)
77 pde.getSolverOptions().setVerbosityOn()
78 pde.setSymmetryOn()
79 #pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
80 psi_ref=pde.getSolution()
81 del pde
82 d_obs=kronecker(DIM)[DIM-1]
83 g_hat=grad(psi_ref)[DIM-1]
84 beta=1/1000000.
85 #
86 # where do we know the gravity:
87 #
88 chi=1.
89 x=Function(domain).getX()
90 dz=H/NE_H
91 H_earth=int(H_earth/dz)*dz # lock to grid
92 chi=whereNegative(abs(x[DIM-1]-(H_earth+dz/2))-dz/2)
93 chi=1
94 #
95 # normalize g_hat (data):
96 #
97 g0=Lsup(chi * g_hat)
98 if not g0 > 0: g0=1.
99 g_hat*=1./g0
100 print "Data normalization factor = %e"%g0
101
102 #===============================================================================================================
103 #
104 # use a linear PDE to solve the problem:
105 #
106 print "====== Use linear PDE ============================================="
107 pde=LinearPDE(domain, numEquations=3, numSolutions=3)
108 A=pde.createCoefficient("A")
109 X=pde.createCoefficient("X")
110 D=pde.createCoefficient("D")
111 q=pde.createCoefficient("q")
112
113
114 A[0,:,0,:]=kronecker(DIM) *beta
115 A[1,:,1,:]=kronecker(DIM)
116 A[2,:,2,:]=kronecker(DIM)
117 A[2,:,1,:]=chi * outer(d_obs, d_obs)
118 D[0,2]=-4*pi/LL**2
119 D[1,0]=-4*pi/LL**2
120 X[2,:]= chi * g_hat * d_obs
121 #q[0]=wherePositive(domain.getX()[DIM-1]-H_earth)
122 #q[0]=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) + whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
123 q[0]=wherePositive(domain.getX()[DIM-1]-H_earth) + whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) + whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
124 #q[0]=whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
125 q[1]=m_psi_ref
126 q[2]=m_psi_ref
127
128 pde.setValue(A=A, D=D, X=X, q=q)
129 pde.getSolverOptions().setVerbosityOn()
130 pde.getSolverOptions().setTolerance(1e-8)
131 pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
132
133 u=pde.getSolution()
134 rho_l, psi_l, =u[0], u[1]
135 print "rho =",rho_l
136 print "psi =",psi_l
137 print "lambda =",u[2]
138
139
140 #=========== This is the same with the variational class ===================================
141 print "====== Use variational problem ============================================="
142 psi_s=Symbol("psi", (), dim=DIM)
143 rho_s=Symbol("rho", (), dim=DIM)
144 #g=Symbol("g", (), dim=DIM)
145
146 v=VariationalProblem(domain, u=psi_s,p=rho_s, debug=VariationalProblem.DEBUG3)
147 v.setValue( H = 0.5*chi*(grad(psi_s)[DIM-1]-g_hat)**2 + 0.5* beta * length(grad(rho_s))**2,
148 X=grad(psi_s), Y=-4*pi*rho_s/LL**2,
149 qp=q[0], q=m_psi_ref)
150 v.getNonlinearPDE().getLinearSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
151
152 rho_v, psi_v, lag=v.getSolution(psi=0, rho=0) # gamma=1 is the interesting case!
153 print "rho =",rho_v
154 print "rho_ref =",rho_ref*g0
155 print "psi =",psi_v
156 print "lambda =",lag
157
158 print "Differences PDE <-> Variational:"
159 print "rho =",Lsup(rho_v-rho_l)/Lsup(rho_l)
160 print "psi =",Lsup(psi_v-psi_l)/Lsup(psi_l)
161 print "Differences to Data :"
162 print "rho =",Lsup(rho_ref-rho_l*g0/LL**2)/Lsup(rho_ref)
163 print "psi =",Lsup(psi_ref-psi_l*g0)/Lsup(psi_ref)
164 print "g =",Lsup(grad(psi_ref-psi_l*g0)[DIM-1])/Lsup(grad(psi_ref)[DIM-1])
165
166
167 #saveVTK("u.vtu", rho_ref=rho_ref, psi_ref=psi_ref, g=grad(psi_ref)[DIM-1], rho=rho_v, psi=psi_v, chi=chi)

  ViewVC Help
Powered by ViewVC 1.1.26