test/python/run_gravity.py


########################################################
#
# Copyright (c) 2003-2012 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

import logging
import os
from esys.escript import unitsSI as U
from esys.escript import inf,sup
from esys.downunder.datasources import SyntheticDataSource,SmoothAnomaly
from esys.downunder.inversions import GravityInversion


try: 
   WORKDIR=os.environ['DOWNUNDER_WORKDIR']
except KeyError:
   WORKDIR='.'

features=[SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km, \
     x=100*U.km, y=3*U.km, depth=25*U.km, rho_inner=200., rho_outer=1e-6),\
          SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km,
     x=400*U.km, y=1*U.km, depth=40*U.km, rho_inner=-200, rho_outer=1e-6)]

logger=logging.getLogger('inv')
logger.setLevel(logging.FATAL)
handler=logging.StreamHandler()
handler.setLevel(logging.FATAL)
logger.addHandler(handler)
source=SyntheticDataSource(DIM=2, NE=20, l=500*U.km, h=60*U.km, features=features)
source.setPadding(5, 0.1)
inv=GravityInversion()
inv.setDataSource(source)
inv.setOutputDirectory(WORKDIR)
inv.setSolverTolerance(1e-5)
inv.setSolverMaxIterations(100)
inv.setSolverOptions(initialHessian=100)
x=source.getDomain().getX()
l0=sup(x[0])-inf(x[0])
l1=sup(x[1])-inf(x[1])
l=max(l0,l1)
G=6.6742e-11
mu=0.5*(l**2*G)**2
inv.setWeights(mu_reg=mu)
rho_new=inv.run()

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=50U.km, ly=20U.km, lz=40*U.km, \
36	x=100U.km, y=3U.km, depth=25*U.km, rho_inner=200., rho_outer=1e-6),\
37	SmoothAnomaly(lx=50U.km, ly=20U.km, lz=40*U.km,
38	x=400U.km, y=1U.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=500U.km, h=60U.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(l2G)**2
59	inv.setWeights(mu_reg=mu)
60	rho_new=inv.run()
61