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

Revision 3877 - (show annotations)
Mon Mar 19 00:52:00 2012 UTC (6 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 5179 byte(s)
```underrelaxtion fixed.
```
 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 x=Function(domain).getX() 89 dz=H/NE_H 90 H_earth=int(H_earth/dz)*dz # lock to grid 91 chi=whereNegative(abs(x[DIM-1]-(H_earth+dz/2))-dz/2) 92 # 93 # normalize g_hat (data): 94 # 95 g0=Lsup(chi * g_hat) 96 if not g0 > 0: g0=1. 97 g_hat*=1./g0 98 print "Data normalization factor = %e"%g0 99 100 #=============================================================================================================== 101 # 102 # use a linear PDE to solve the problem: 103 # 104 print "====== Use linear PDE =============================================" 105 pde=LinearPDE(domain, numEquations=3, numSolutions=3) 106 A=pde.createCoefficient("A") 107 X=pde.createCoefficient("X") 108 D=pde.createCoefficient("D") 109 q=pde.createCoefficient("q") 110 111 112 A[0,:,0,:]=kronecker(DIM) *beta 113 A[1,:,1,:]=kronecker(DIM) 114 A[2,:,2,:]=kronecker(DIM) 115 A[2,:,1,:]=chi * outer(d_obs, d_obs) 116 D[0,2]=-4*pi/LL**2 117 D[1,0]=-4*pi/LL**2 118 X[2,:]= chi * g_hat * d_obs 119 #q[0]=wherePositive(domain.getX()[DIM-1]-H_earth) 120 #q[0]=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) + whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1])) 121 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])) 122 #q[0]=whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1])) 123 q[1]=m_psi_ref 124 q[2]=m_psi_ref 125 126 pde.setValue(A=A, D=D, X=X, q=q) 127 pde.getSolverOptions().setVerbosityOn() 128 pde.getSolverOptions().setTolerance(1e-8) 129 pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT) 130 131 u=pde.getSolution() 132 rho_l, psi_l, =u[0], u[1] 133 print "rho =",rho_l 134 print "psi =",psi_l 135 print "lambda =",u[2] 136 137 138 #=========== This is the same with the variational class =================================== 139 print "====== Use variational problem =============================================" 140 psi_s=Symbol("psi", (), dim=DIM) 141 rho_s=Symbol("rho", (), dim=DIM) 142 #g=Symbol("g", (), dim=DIM) 143 144 v=VariationalProblem(domain, u=psi_s,p=rho_s, debug=VariationalProblem.DEBUG3) 145 v.setValue( H = 0.5*chi*(grad(psi_s)[DIM-1]-g_hat)**2 + 0.5* beta * length(grad(rho_s))**2, 146 X=grad(psi_s), Y=-4*pi*rho_s/LL**2, 147 qp=q[0], q=m_psi_ref) 148 v.getNonlinearPDE().getLinearSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT) 149 150 rho_v, psi_v, lag=v.getSolution(psi=0, rho=0) # gamma=1 is the interesting case! 151 print "rho =",rho_v 152 print "rho_ref =",rho_ref*g0 153 print "psi =",psi_v 154 print "lambda =",lag 155 156 print "Differences PDE <-> Variational:" 157 print "rho =",Lsup(rho_v-rho_l)/Lsup(rho_l) 158 print "psi =",Lsup(psi_v-psi_l)/Lsup(psi_l) 159 print "Differences to Data :" 160 print "rho =",Lsup(rho_ref-rho_l*g0/LL**2)/Lsup(rho_ref) 161 print "psi =",Lsup(psi_ref-psi_l*g0)/Lsup(psi_ref) 162 print "g =",Lsup(chi*grad(psi_ref-psi_l*g0)[DIM-1])/Lsup(chi*grad(psi_ref)[DIM-1]) 163 164 165 #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)