/[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 4267 - (show annotations)
Fri Mar 1 04:50:54 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 4143 byte(s)
Updated inversion example scripts, added ER Mapper data, checked results
and corrected some issues. Cookbook does not reflect the changes yet.

1 ##############################################################################
2 #
3 # Copyright (c) 2012-2013 by 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 since 2012 by School of Earth Sciences
12 #
13 ##############################################################################
14
15 """
16 Advanced 3D gravity/magnetic joint inversion example without using any
17 inversion drivers
18 """
19
20 __copyright__="""Copyright (c) 2012-2013 by University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Open Software License version 3.0
24 http://www.opensource.org/licenses/osl-3.0.php"""
25 __url__="https://launchpad.net/escript-finley"
26
27 # Import required modules
28 import numpy as np
29 from esys.downunder import *
30 from esys.escript import unitsSI as U
31 from esys.escript import *
32 from esys.weipa import *
33
34 # Set parameters
35 MAGNETIC_DATASET = 'data/MagneticSmall.nc'
36 MAG_UNITS = U.Nano * U.V * U.sec / (U.m**2)
37 GRAVITY_DATASET = 'data/GravitySmall.nc'
38 GRAV_UNITS = 1e-6 * U.m/(U.sec**2)
39 latitude=-20.5
40
41 thickness = 40. * U.km # below surface
42 l_air = 6. * U.km # above surface
43 n_cells_v = 25 # number of cells in vertical direction
44
45 # apply 20% padding
46 PAD_X = 0.2
47 PAD_Y = 0.2
48
49 MU_GRAVITY = 10.
50 MU_MAGNETIC = 0.1
51
52 # read data:
53 source_g=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS)
54 source_m=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS)
55 B_b=simpleGeoMagneticFluxDensity(latitude=latitude)
56
57 # create domain:
58 db=DomainBuilder(dim=3)
59 db.addSource(source_g)
60 db.addSource(source_m)
61 db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)
62 db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)
63 db.fixDensityBelow(depth=thickness)
64 db.fixSusceptibilityBelow(depth=thickness)
65
66 dom=db.getDomain()
67 DIM=dom.getDim()
68
69 # create mappings with standard parameters
70 rho_mapping=DensityMapping(dom)
71 k_mapping=SusceptibilityMapping(dom)
72
73 # create regularization with two level set functions:
74 reg_mask=Data(0.,(2,), Solution(dom))
75 reg_mask[0] = db.getSetDensityMask() # mask for locations where m[0]~rho is known
76 reg_mask[1] = db.getSetSusceptibilityMask() # mask for locations where m[0]~k is known
77 regularization=Regularization(dom, numLevelSets=2,
78 w1=np.ones((2,DIM)), # consider gradient terms
79 wc=[[0,1],[0,0]], # and cross-gradient term
80 location_of_set_m=reg_mask)
81
82 # create forward model for gravity
83 # get data with deviation
84 g,sigma_g=db.getGravitySurveys()[0]
85 # turn the scalars into vectors (vertical direction)
86 d=kronecker(DIM)[DIM-1]
87 w=safeDiv(1., sigma_g)
88
89 gravity_model=GravityModel(dom, w*d, g*d)
90 gravity_model.rescaleWeights(rho_scale=rho_mapping.getTypicalDerivative())
91
92 # create forward model for magnetic
93 d=normalize(np.array(B_b)) # direction of measurement
94 B,sigma_B=db.getMagneticSurveys()[0]
95 w=safeDiv(1., sigma_B)
96
97 magnetic_model=MagneticModel(dom, w*d, B*d, B_b)
98 magnetic_model.rescaleWeights(k_scale=k_mapping.getTypicalDerivative())
99
100
101 # finally we can set up the cost function:
102 cf=InversionCostFunction(regularization,
103 ((rho_mapping,0), (k_mapping, 1)),
104 ((gravity_model,0), (magnetic_model,1)) )
105
106 cf.setTradeOffFactorsModels([MU_GRAVITY, MU_MAGNETIC])
107
108 # sun solver:
109 solver=MinimizerLBFGS()
110 solver.setCostFunction(cf)
111 solver.setTolerance(1e-4)
112 solver.setMaxIterations(50)
113 solver.run(Data(0.,(2,),Solution(dom)))
114 m=solver.getResult()
115 density, susceptibility = cf.getProperties(m)
116
117
118 # write everything to file:
119 saveSilo("result_gravmag.silo",
120 density=density, susceptability=susceptibility,
121 g_data=g, sigma_g=sigma_g, B_data=B, sigma_B=sigma_B)
122 saveVTK("result_gravmag.vtu",
123 density=density, susceptability=susceptibility,
124 g_data=g, sigma_g=sigma_g, B_data=B, sigma_B=sigma_B)
125
126 print("All done. Have a nice day!")
127

  ViewVC Help
Powered by ViewVC 1.1.26