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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.4266  
changed lines
  Added in v.4267

  ViewVC Help
Powered by ViewVC 1.1.26