/[escript]/trunk/doc/examples/inversion/gravmag_nodriver.py
ViewVC logotype

Contents of /trunk/doc/examples/inversion/gravmag_nodriver.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 months ago) by sshaw
File MIME type: text/x-python
File size: 5438 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2012-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 """
18 Advanced 3D gravity/magnetic joint inversion example without using any
19 inversion drivers
20 """
21
22 __copyright__="""Copyright (c) 2012-2015 by The University of Queensland
23 http://www.uq.edu.au
24 Primary Business: Queensland, Australia"""
25 __license__="""Licensed under the Open Software License version 3.0
26 http://www.opensource.org/licenses/osl-3.0.php"""
27 __url__="https://launchpad.net/escript-finley"
28
29 # Import required modules
30 import numpy as np
31 from esys.downunder import *
32 from esys.escript import unitsSI as U
33 from esys.escript import *
34 from esys.weipa import *
35
36 haveNetcdf=(getEscriptParamInt("NETCDF_BUILD",0)==1)
37 try:
38 import pyproj
39 havePyProj=True
40 except ImportError:
41 havePyProj=False
42
43 # Set parameters
44 MAGNETIC_DATASET = 'data/MagneticSmall.nc'
45 MAG_UNITS = U.Nano * U.Tesla
46 GRAVITY_DATASET = 'data/GravitySmall.nc'
47 GRAV_UNITS = 1e-6 * U.m/(U.sec**2)
48 # background magnetic field components (B_East, B_North, B_Vertical)
49 B_b = [2201.*U.Nano*U.Tesla, 31232.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]
50
51 thickness = 40. * U.km # below surface
52 l_air = 6. * U.km # above surface
53 n_cells_v = 25 # number of cells in vertical direction
54
55 # apply 20% padding
56 PAD_X = 0.2
57 PAD_Y = 0.2
58
59 MU_GRAVITY = 10.
60 MU_MAGNETIC = 0.1
61
62 COORDINATES=CartesianReferenceSystem()
63
64 def work():
65 # read data:
66 source_g=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS, reference_system=COORDINATES)
67 source_m=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS, reference_system=COORDINATES)
68
69 # create domain:
70 db=DomainBuilder(dim=3, reference_system=COORDINATES)
71 db.addSource(source_g)
72 db.addSource(source_m)
73 db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)
74 db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)
75 db.fixDensityBelow(depth=thickness)
76 db.fixSusceptibilityBelow(depth=thickness)
77
78 dom=db.getDomain()
79 DIM=dom.getDim()
80
81 # create mappings with standard parameters
82 rho_mapping=DensityMapping(dom)
83 k_mapping=SusceptibilityMapping(dom)
84
85 # create regularization with two level set functions:
86 reg_mask=Data(0.,(2,), Solution(dom))
87 reg_mask[0] = db.getSetDensityMask() # mask for locations where m[0]~rho is known
88 reg_mask[1] = db.getSetSusceptibilityMask() # mask for locations where m[0]~k is known
89 regularization=Regularization(dom, numLevelSets=2,
90 w1=np.ones((2,DIM)), # consider gradient terms
91 wc=[[0,1],[0,0]], # and cross-gradient term
92 coordinates=COORDINATES,
93 location_of_set_m=reg_mask)
94
95 # create forward model for gravity
96 # get data with deviation
97 g,sigma_g=db.getGravitySurveys()[0]
98 # turn the scalars into vectors (vertical direction)
99 d=kronecker(DIM)[DIM-1] # == (0 0 1)
100 w=safeDiv(1., sigma_g)
101
102 # multipling by d extracts only the z component (since we only measure in vertical)
103 gravity_model=GravityModel(dom, w*d, g*d, coordinates=COORDINATES)
104 gravity_model.rescaleWeights(rho_scale=rho_mapping.getTypicalDerivative())
105
106 # create forward model for magnetic
107 d=normalize(np.array(B_b)) # direction of measurement
108 B,sigma_B=db.getMagneticSurveys()[0]
109 w=safeDiv(1., sigma_B)
110
111 magnetic_model=MagneticModel(dom, w*d, B*d, B_b, coordinates=COORDINATES)
112 # or:
113 # magnetic_model=SelfDemagnetizationModel(dom, w*d, B*d, B_b, coordinates=COORDINATES)
114 magnetic_model.rescaleWeights(k_scale=k_mapping.getTypicalDerivative())
115
116
117 # finally we can set up the cost function:
118 cf=InversionCostFunction(regularization,
119 ((rho_mapping,0), (k_mapping, 1)),
120 ((gravity_model,0), (magnetic_model,1)) )
121
122 cf.setTradeOffFactorsModels([MU_GRAVITY, MU_MAGNETIC])
123
124 # sun solver:
125 solver=MinimizerLBFGS()
126 solver.setCostFunction(cf)
127 solver.setTolerance(1e-4)
128 solver.setMaxIterations(50)
129 solver.run(Data(0.,(2,),Solution(dom)))
130 m=solver.getResult()
131 density, susceptibility = cf.getProperties(m)
132
133
134 # write everything to file:
135 saveSilo("result_gravmag.silo",
136 density=density, susceptability=susceptibility,
137 g_data=g, sigma_g=sigma_g, B_data=B, sigma_B=sigma_B)
138 saveVTK("result_gravmag.vtu",
139 density=density, susceptability=susceptibility,
140 g_data=g, sigma_g=sigma_g, B_data=B, sigma_B=sigma_B)
141
142 print("All done. Have a nice day!")
143
144 try:
145 import esys.ripley
146 HAVE_RIPLEY = True
147 except ImportError:
148 HAVE_RIPLEY = False
149
150 if 'NetCdfData' not in dir():
151 print("This example requires scipy's netcdf support which does not appear to be installed.")
152 elif not HAVE_RIPLEY:
153 print("Ripley module not available")
154 elif not haveNetcdf:
155 print("netCDF not available.")
156 elif not havePyProj:
157 print("This example requires pyproj.")
158 else:
159 work()
160

  ViewVC Help
Powered by ViewVC 1.1.26