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

  ViewVC Help
Powered by ViewVC 1.1.26