| 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=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= m_psi + whereZero(x_cord[i]-inf(x_cord[i])) + whereZero(x_cord[i]-sup(x_cord[i])) |
| 65 |
#m_rho=wherePositive(domain.getX()[DIM-1]-H_earth) |
| 66 |
#m_rho=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) + whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1])) |
| 67 |
m_rho=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])) |
| 68 |
#m_rho=whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1])) |
| 69 |
|
| 70 |
|
| 71 |
# create test density: |
| 72 |
rho_ref= 0 |
| 73 |
for f in feastures: |
| 74 |
m=f.getMask(x_cord) |
| 75 |
rho_ref = rho_ref * (1-m) + f.rho * m |
| 76 |
|
| 77 |
#rho_ref=sup(x_cord[DIM-1])-x_cord[DIM-1] for testing! |
| 78 |
|
| 79 |
# get the reference potential: |
| 80 |
pde=LinearSinglePDE(domain) |
| 81 |
pde.setValue(A=kronecker(domain), Y=4*pi*rho_ref, q=m_psi) |
| 82 |
pde.getSolverOptions().setVerbosityOn() |
| 83 |
pde.setSymmetryOn() |
| 84 |
#pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT) |
| 85 |
psi_ref=pde.getSolution() |
| 86 |
del pde |
| 87 |
d_obs=kronecker(DIM)[DIM-1] |
| 88 |
g_hat=grad(psi_ref)[DIM-1] |
| 89 |
beta=1/1000000. |
| 90 |
#beta=1/100. |
| 91 |
# |
| 92 |
# where do we know the gravity: |
| 93 |
# |
| 94 |
chi=1. |
| 95 |
x=Function(domain).getX() |
| 96 |
dz=H/NE_H |
| 97 |
H_earth=int(H_earth/dz)*dz # lock to grid |
| 98 |
chi=whereNegative(abs(x[DIM-1]-(H_earth+dz/2))-dz/2) |
| 99 |
# |
| 100 |
# normalize g_hat (data): |
| 101 |
# |
| 102 |
g0=Lsup(chi * g_hat) |
| 103 |
if not g0 > 0: g0=1. |
| 104 |
g_hat*=1./g0 |
| 105 |
print "Data normalization factor = %e"%g0 |
| 106 |
|
| 107 |
#=========== This is the same with the variational class =================================== |
| 108 |
print "====== Use variational problem =============================================" |
| 109 |
psi_s=Symbol("psi", (), dim=DIM) |
| 110 |
rho_s=Symbol("rho", (), dim=DIM) |
| 111 |
gamma_s=Symbol("gamma", (), dim=DIM) |
| 112 |
gamma_s=0.5 |
| 113 |
gamma_s=0.6 |
| 114 |
#g=Symbol("g", (), dim=DIM) |
| 115 |
|
| 116 |
v=VariationalProblem(domain, u=psi_s,p=rho_s, debug=VariationalProblem.DEBUG3) |
| 117 |
v.setValue( H = 0.5*chi*(grad(psi_s)[DIM-1]-g_hat)**2 + beta/gamma_s * (length(grad(rho_s))**2 + EPSILON**2)**gamma_s, |
| 118 |
X=grad(psi_s), Y=-4*pi*rho_s/LL**2, |
| 119 |
qp=m_rho, q=m_psi) |
| 120 |
v.getNonlinearPDE().getLinearSolverOptions().setSolverMethod(v.getNonlinearPDE().getLinearSolverOptions().DIRECT) |
| 121 |
|
| 122 |
rho_v, psi_v, lag=v.getSolution(psi=0, rho=1, gamma=2) # gamma=0.5 is the interesting case! |
| 123 |
print "rho =",rho_v |
| 124 |
print "rho_ref =",rho_ref*g0 |
| 125 |
print "psi =",psi_v |
| 126 |
print "lambda =",lag |
| 127 |
|
| 128 |
print "Differences to Data :" |
| 129 |
print "rho =",Lsup(rho_ref-rho_v*g0/LL**2)/Lsup(rho_ref) |
| 130 |
print "psi =",Lsup(psi_ref-psi_v*g0)/Lsup(psi_ref) |
| 131 |
print "g =",Lsup(chi*grad(psi_ref-psi_v*g0)[DIM-1])/Lsup(chi*grad(psi_ref)[DIM-1]) |
| 132 |
|
| 133 |
|
| 134 |
|
| 135 |
#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) |
Parent Directory
| 
Revision Log