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

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

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

revision 4266 by caltinay, Wed Feb 27 00:26:12 2013 UTC revision 4267 by caltinay, Fri Mar 1 04:50:54 2013 UTC
# Line 24  __url__="https://launchpad.net/escript-f Line 24  __url__="https://launchpad.net/escript-f
24    
25  # Import required modules  # Import required modules
26  from esys.downunder import *  from esys.downunder import *
 from esys.downunder.minimizers import MinimizerMaxIterReached  
27  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
28  from esys.weipa import saveSilo  from esys.escript import saveDataCSV
29    from esys.weipa import *
30    
31  # Set parameters  # Set parameters
32  MAGNETIC_DATASET = 'magnetic_anomaly.nc'  MAGNETIC_DATASET = 'data/MagneticSmall.nc'
33  GRAVITY_DATASET = 'gravity_anomaly.nc'  MAG_UNITS = U.Nano * U.V * U.sec / (U.m**2)
34  latitude = -28.5  GRAVITY_DATASET = 'data/GravitySmall.nc'
35    GRAV_UNITS = 1e-6 * U.m/(U.sec**2)
36    latitude = -20.5
37  PAD_X = 0.2  PAD_X = 0.2
38  PAD_Y = 0.2  PAD_Y = 0.2
39  thickness = 40. * U.km  thickness = 40. * U.km
# Line 42  mu_magnetic = 0.1 Line 44  mu_magnetic = 0.1
44    
45  # Setup and run the inversion  # Setup and run the inversion
46  B_b=simpleGeoMagneticFluxDensity(latitude=latitude)  B_b=simpleGeoMagneticFluxDensity(latitude=latitude)
47  grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET)  grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS)
48  mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET)  mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS)
49  db=DomainBuilder(dim=3)  db=DomainBuilder(dim=3)
50  db.addSource(grav_source)  db.addSource(grav_source)
51  db.addSource(mag_source)  db.addSource(mag_source)
# Line 60  inv.setup(db) Line 62  inv.setup(db)
62  inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic])  inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic])
63  inv.getCostFunction().setTradeOffFactorsRegularization(mu = [1.,1.], mu_c=1.)  inv.getCostFunction().setTradeOffFactorsRegularization(mu = [1.,1.], mu_c=1.)
64    
65  print("Starting inversion, please stand by...")  density, susceptibility = inv.run()
 try:  
     density, susceptibility = inv.run()  
 except MinimizerMaxIterReached as e:  
     print(e)  
     density, susceptibility = inv.p  
   
66  print("density = %s"%density)  print("density = %s"%density)
67  print("susceptibility = %s"%susceptibility)  print("susceptibility = %s"%susceptibility)
68    
69  # Save results  g, wg = db.getGravitySurveys()[0]
70  gravSurveys = db.getGravitySurveys()  B, wB = db.getMagneticSurveys()[0]
71  magSurveys = db.getMagneticSurveys()  saveSilo("result_gravmag.silo", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B, magnetic_weight=wB)
72  g, wg =  gravSurveys[0]  print("Results saved in result_gravmag.silo")
73  B, wB =  magSurveys[0]  
74  saveSilo("result_gravmag.silo", density=density, gravity_anomaly=g[2], gravity_weight=wg[2], susceptibility=susceptibility, magnetic_anomaly=B[2], magnetic_weight=wB[2])  saveVTK("result_gravmag.vtu", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B, magnetic_weight=wB)
75  print("Results saved in result_gravmag.silo. Good bye.")  print("Results saved in result_gravmag.vtu")
76    
77    saveDataCSV("result_gravmag.csv", density=density, susceptibility=susceptibility, x=susceptibility.getFunctionSpace().getX())
78    print("Results saved in result_gravmag.csv")
79    
80    print("All done. Have a nice day!")
81    

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

  ViewVC Help
Powered by ViewVC 1.1.26