test/python/run_magnetic.py


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

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
http://www.uq.edu.au
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.downunder import *
from esys.escript import unitsSI as U
from esys.weipa import saveSilo

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

features=[SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \
     x=8*U.km, y=3*U.km, depth=2.5*U.km, v_inner=2., v_outer=1e-6),\
          SmoothAnomaly(lx=25*U.km, ly=20*U.km, lz=20*U.km,
     x=30*U.km, y=1*U.km, depth=18*U.km, v_inner=-1., v_outer=1e-6),\
          SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \
     x=68*U.km, y=3*U.km, depth=5*U.km, v_inner=20., v_outer=1e-6)]


logger=logging.getLogger('inv')
logger.setLevel(logging.INFO)
handler=logging.StreamHandler()
handler.setLevel(logging.INFO)
logger.addHandler(handler)
source=SyntheticData(DataSource.MAGNETIC, DIM=2, NE=30, l=100*U.km, features=features)
domainbuilder=DomainBuilder(dim=2)
domainbuilder.addSource(source)
domainbuilder.setPadding(10)
domainbuilder.setVerticalExtents(depth=30*U.km, air_layer=10*U.km, num_cells=16)

inv=MagneticInversion()
inv.setSolverTolerance(1e-4)
inv.setSolverMaxIterations(10)
inv.setSolverOptions(initialHessian=100)
#inv.setWeights(mu_reg=1e-4)
inv.setup(domainbuilder)

k_new=inv.run()
B, chi = inv.getForwardModel().getSurvey(0)
saveSilo(os.path.join(WORKDIR, 'maginv'), sus=k_new, sus_mask=inv.getRegularization().location_of_set_m, sus_ref=source.getReferenceSusceptibility(), B=B, chi=chi)

1
2	##############################################################################
3	#
4	# Copyright (c) 2003-2012 by University of Queensland
5	# http://www.uq.edu.au
6	#
7	# Primary Business: Queensland, Australia
8	# Licensed under the Open Software License version 3.0
9	# http://www.opensource.org/licenses/osl-3.0.php
10	#
11	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	# Development since 2012 by School of Earth Sciences
13	#
14	##############################################################################
15
16	__copyright__="""Copyright (c) 2003-2012 by University of Queensland
17	http://www.uq.edu.au
18	Primary Business: Queensland, Australia"""
19	__license__="""Licensed under the Open Software License version 3.0
20	http://www.opensource.org/licenses/osl-3.0.php"""
21	__url__="https://launchpad.net/escript-finley"
22
23	import logging
24	import os
25	from esys.downunder import *
26	from esys.escript import unitsSI as U
27	from esys.weipa import saveSilo
28
29	try:
30	WORKDIR=os.environ['DOWNUNDER_WORKDIR']
31	except KeyError:
32	WORKDIR='.'
33
34	features=[SmoothAnomaly(lx=30U.km, ly=20U.km, lz=18.*U.km, \
35	x=8U.km, y=3U.km, depth=2.5*U.km, v_inner=2., v_outer=1e-6),\
36	SmoothAnomaly(lx=25U.km, ly=20U.km, lz=20*U.km,
37	x=30U.km, y=1U.km, depth=18*U.km, v_inner=-1., v_outer=1e-6),\
38	SmoothAnomaly(lx=30U.km, ly=20U.km, lz=18.*U.km, \
39	x=68U.km, y=3U.km, depth=5*U.km, v_inner=20., v_outer=1e-6)]
40
41
42	logger=logging.getLogger('inv')
43	logger.setLevel(logging.INFO)
44	handler=logging.StreamHandler()
45	handler.setLevel(logging.INFO)
46	logger.addHandler(handler)
47	source=SyntheticData(DataSource.MAGNETIC, DIM=2, NE=30, l=100*U.km, features=features)
48	domainbuilder=DomainBuilder(dim=2)
49	domainbuilder.addSource(source)
50	domainbuilder.setPadding(10)
51	domainbuilder.setVerticalExtents(depth=30U.km, air_layer=10U.km, num_cells=16)
52
53	inv=MagneticInversion()
54	inv.setSolverTolerance(1e-4)
55	inv.setSolverMaxIterations(10)
56	inv.setSolverOptions(initialHessian=100)
57	#inv.setWeights(mu_reg=1e-4)
58	inv.setup(domainbuilder)
59
60	k_new=inv.run()
61	B, chi = inv.getForwardModel().getSurvey(0)
62	saveSilo(os.path.join(WORKDIR, 'maginv'), sus=k_new, sus_mask=inv.getRegularization().location_of_set_m, sus_ref=source.getReferenceSusceptibility(), B=B, chi=chi)
63