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

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