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

Diff of /branches/doubleplusgood/doc/examples/inversion/gravmag_nodriver.py

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

revision 4344 by jfenwick, Wed Feb 27 03:42:40 2013 UTC revision 4345 by jfenwick, Fri Mar 29 07:09:41 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    # background magnetic field components (B_North, B_East, B_Vertical)
40    B_b = [31232.*U.Nano*U.Tesla, 2201.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]
41    
42  thickness = 40. * U.km # below surface  thickness = 40. * U.km # below surface
43  l_air = 6. * U.km      # above surface  l_air = 6. * U.km      # above surface
44  n_cells_v = 25         # number of cells in vertical direction    n_cells_v = 25         # number of cells in vertical direction
45    
46  # apply 20% padding  # apply 20% padding
47  PAD_X = 0.2  PAD_X = 0.2
48  PAD_Y = 0.2  PAD_Y = 0.2
49    
50  mu = 10.  MU_GRAVITY = 10.
51    MU_MAGNETIC = 0.1
52    
53  # read data:  def work():
54  source_g=NetCdfData(NetCdfData.GRAVITY, DATASET_GRAV)    # read data:
55  source_m=NetCdfData(NetCdfData.MAGNETIC, DATASET_MAG)    source_g=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS)
56  B_b=simpleGeoMagneticFluxDensity(latitude=latitude)    source_m=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS)
57    
58      # create domain:
59  # create domain:    db=DomainBuilder(dim=3)
60  db=DomainBuilder(dim=3)    db.addSource(source_g)
61  db.addSource(source_g)    db.addSource(source_m)
62  db.addSource(source_m)    db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)
63  db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)    db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)
64  db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)    db.fixDensityBelow(depth=thickness)
65  db.fixDensityBelow(depth=thickness)    db.fixSusceptibilityBelow(depth=thickness)
66    
67  dom=db.getDomain()    dom=db.getDomain()
68  DIM=dom.getDim()    DIM=dom.getDim()
69    
70  # create mappings with standard parameters    # create mappings with standard parameters
71  rho_mapping=DensityMapping(dom)    rho_mapping=DensityMapping(dom)
72  k_mapping=SusceptibilityMapping(dom)    k_mapping=SusceptibilityMapping(dom)
73            
74  # create regularization with two level set functions:    # create regularization with two level set functions:
75  reg_mask=Data(0.,(2,), Solution(dom))      reg_mask=Data(0.,(2,), Solution(dom))
76  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
77  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
78  regularization=Regularization(dom, numLevelSets=2,    regularization=Regularization(dom, numLevelSets=2,
79                                 w1=np.ones((2,DIM)),      # consider gradient terms                                 w1=np.ones((2,DIM)), # consider gradient terms
80                                 wc=[[0,0],[1,0]],         # and cross gradient term                                 wc=[[0,1],[0,0]],    # and cross-gradient term
81                                 location_of_set_m=reg_mask)                                 location_of_set_m=reg_mask)
                                 
 # create forward model for gravity  
 # get data with deviation  
 g,sigma_g=domainbuilder.getGravitySurveys()[0]  
 # turn the scalars into vectors (vertical direction)  
 d=kronecker(DIM)[DIM-1]  
 w=safeDiv(1., sigma_g)  
   
 gravity_model=GravityModel(dom, w*d, g*d)  
 gravity_model.rescaleWeights(rho_scale=rho_mapping.getTypicalDerivative())  
   
 # create forward model for magnetic  
 d=normalize(B_b) # direction of measurement  
 B,sigma_B=domainbuilder.getMagneticSurveys()[0]  
 w=safeDiv(1., sigma_B)  
82    
83  magnetic_model=MagneticModel(dom, w*d, B*d, B_b)    # create forward model for gravity
84  magnetic_model.rescaleWeights(k_scale=k_mapping.getTypicalDerivative())    # get data with deviation
85      g,sigma_g=db.getGravitySurveys()[0]
86      # turn the scalars into vectors (vertical direction)
87      d=kronecker(DIM)[DIM-1]
88      w=safeDiv(1., sigma_g)
89    
90      gravity_model=GravityModel(dom, w*d, g*d)
91      gravity_model.rescaleWeights(rho_scale=rho_mapping.getTypicalDerivative())
92    
93      # create forward model for magnetic
94      d=normalize(np.array(B_b)) # direction of measurement
95      B,sigma_B=db.getMagneticSurveys()[0]
96      w=safeDiv(1., sigma_B)
97    
98      magnetic_model=MagneticModel(dom, w*d, B*d, B_b)
99      magnetic_model.rescaleWeights(k_scale=k_mapping.getTypicalDerivative())
100    
101    
102  # finally we can set up the cost function:    # finally we can set up the cost function:
103  cf=InversionCostFunction(regularization,    cf=InversionCostFunction(regularization,
104               ((rho_mapping,0), (k_mapping, 1)),               ((rho_mapping,0), (k_mapping, 1)),
105               ((gravity_model,0), (magnetic_model,1)) )               ((gravity_model,0), (magnetic_model,1)) )
106    
107  cf.setTradeOffFactorsModels(mu)    cf.setTradeOffFactorsModels([MU_GRAVITY, MU_MAGNETIC])
108    
109  # sun solver:    # sun solver:
110  solver=MinimizerLBFGS()    solver=MinimizerLBFGS()
111  solver.setSolverTolerance(1e-4)    solver.setCostFunction(cf)
112  solver.setSolverMaxIterations(50)    solver.setTolerance(1e-4)
113  m=solver.run(Data(0.,(2,),Solution(dom)))    solver.setMaxIterations(50)
114  density, susceptibility =getProperties(m)    solver.run(Data(0.,(2,),Solution(dom)))
115      m=solver.getResult()
116      density, susceptibility = cf.getProperties(m)
117    
118    
119  # write everything to file:    # write everything to file:
120  saveSilo("results.silo",    saveSilo("result_gravmag.silo",
121           density=density, susceptability=susceptibility,           density=density, susceptability=susceptibility,
122           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)
123      saveVTK("result_gravmag.vtu",
124             density=density, susceptability=susceptibility,
125             g_data=g, sigma_g=sigma_g, B_data=B, sigma_B=sigma_B)
126    
127      print("All done. Have a nice day!")
128    
129    if 'NetCdfData' in dir():
130      work()
131    else:
132      print("This example requires scipy's netcdf support which does not appear to be installed.")
133    

Legend:
Removed from v.4344  
changed lines
  Added in v.4345

  ViewVC Help
Powered by ViewVC 1.1.26