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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26