1 |
|
2 |
######################################################## |
3 |
# |
4 |
# Copyright (c) 2003-2012 by University of Queensland |
5 |
# Earth Systems Science Computational Center (ESSCC) |
6 |
# http://www.uq.edu.au/esscc |
7 |
# |
8 |
# Primary Business: Queensland, Australia |
9 |
# Licensed under the Open Software License version 3.0 |
10 |
# http://www.opensource.org/licenses/osl-3.0.php |
11 |
# |
12 |
######################################################## |
13 |
|
14 |
__copyright__="""Copyright (c) 2003-2012 by University of Queensland |
15 |
Earth Systems Science Computational Center (ESSCC) |
16 |
http://www.uq.edu.au/esscc |
17 |
Primary Business: Queensland, Australia""" |
18 |
__license__="""Licensed under the Open Software License version 3.0 |
19 |
http://www.opensource.org/licenses/osl-3.0.php""" |
20 |
__url__="https://launchpad.net/escript-finley" |
21 |
|
22 |
import logging |
23 |
import os |
24 |
from esys.escript import unitsSI as U |
25 |
from esys.escript import inf,sup |
26 |
from esys.downunder.datasources import SyntheticDataSource,SmoothAnomaly |
27 |
from esys.downunder.inversions import GravityInversion |
28 |
|
29 |
|
30 |
try: |
31 |
WORKDIR=os.environ['DOWNUNDER_WORKDIR'] |
32 |
except KeyError: |
33 |
WORKDIR='.' |
34 |
|
35 |
features=[SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km, \ |
36 |
x=100*U.km, y=3*U.km, depth=25*U.km, rho_inner=200., rho_outer=1e-6),\ |
37 |
SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km, |
38 |
x=400*U.km, y=1*U.km, depth=40*U.km, rho_inner=-200, rho_outer=1e-6)] |
39 |
|
40 |
logger=logging.getLogger('inv') |
41 |
logger.setLevel(logging.FATAL) |
42 |
handler=logging.StreamHandler() |
43 |
handler.setLevel(logging.FATAL) |
44 |
logger.addHandler(handler) |
45 |
source=SyntheticDataSource(DIM=2, NE=20, l=500*U.km, h=60*U.km, features=features) |
46 |
source.setPadding(5, 0.1) |
47 |
inv=GravityInversion() |
48 |
inv.setDataSource(source) |
49 |
inv.setOutputDirectory(WORKDIR) |
50 |
inv.setSolverTolerance(1e-5) |
51 |
inv.setSolverMaxIterations(100) |
52 |
inv.setSolverOptions(initialHessian=100) |
53 |
x=source.getDomain().getX() |
54 |
l0=sup(x[0])-inf(x[0]) |
55 |
l1=sup(x[1])-inf(x[1]) |
56 |
l=max(l0,l1) |
57 |
G=6.6742e-11 |
58 |
mu=0.5*(l**2*G)**2 |
59 |
inv.setWeights(mu_reg=mu) |
60 |
rho_new=inv.run() |
61 |
|