/[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 4309 by caltinay, Thu Mar 7 05:26:44 2013 UTC revision 4310 by jfenwick, Tue Mar 12 03:23:26 2013 UTC
# Line 43  n_cells_v = 25 Line 43  n_cells_v = 25
43  mu_gravity = 10.  mu_gravity = 10.
44  mu_magnetic = 0.1  mu_magnetic = 0.1
45    
46  # Setup and run the inversion  def work():
47  grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS)    # Setup and run the inversion
48  mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS)    grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS)
49  db=DomainBuilder(dim=3)    mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS)
50  db.addSource(grav_source)    db=DomainBuilder(dim=3)
51  db.addSource(mag_source)    db.addSource(grav_source)
52  db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)    db.addSource(mag_source)
53  db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)    db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)
54  db.setBackgroundMagneticFluxDensity(B_b)    db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)
55  db.fixDensityBelow(depth=thickness)    db.setBackgroundMagneticFluxDensity(B_b)
56  db.fixSusceptibilityBelow(depth=thickness)    db.fixDensityBelow(depth=thickness)
57      db.fixSusceptibilityBelow(depth=thickness)
58  inv=JointGravityMagneticInversion()  
59  inv.setSolverTolerance(1e-4)    inv=JointGravityMagneticInversion()
60  inv.setSolverMaxIterations(50)    inv.setSolverTolerance(1e-4)
61  inv.setup(db)    inv.setSolverMaxIterations(50)
62  inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic])    inv.setup(db)
63  inv.getCostFunction().setTradeOffFactorsRegularization(mu = [1.,1.], mu_c=1.)    inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic])
64      inv.getCostFunction().setTradeOffFactorsRegularization(mu = [1.,1.], mu_c=1.)
65  density, susceptibility = inv.run()  
66  print("density = %s"%density)    density, susceptibility = inv.run()
67  print("susceptibility = %s"%susceptibility)    print("density = %s"%density)
68      print("susceptibility = %s"%susceptibility)
69  g, wg = db.getGravitySurveys()[0]  
70  B, wB = db.getMagneticSurveys()[0]    g, wg = db.getGravitySurveys()[0]
71  saveSilo("result_gravmag.silo", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B, magnetic_weight=wB)    B, wB = db.getMagneticSurveys()[0]
72  print("Results saved in result_gravmag.silo")    saveSilo("result_gravmag.silo", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B,   magnetic_weight=wB)
73      print("Results saved in result_gravmag.silo")
74  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.vtu")    saveVTK("result_gravmag.vtu", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B,   magnetic_weight=wB)
76      print("Results saved in result_gravmag.vtu")
77  saveDataCSV("result_gravmag.csv", density=density, susceptibility=susceptibility, x=susceptibility.getFunctionSpace().getX())  
78  print("Results saved in result_gravmag.csv")    saveDataCSV("result_gravmag.csv", density=density, susceptibility=susceptibility, x=susceptibility.getFunctionSpace().getX())
79      print("Results saved in result_gravmag.csv")
80  print("All done. Have a nice day!")  
81      print("All done. Have a nice day!")
82    
83    if 'NetCdfData' in dir():
84      work()
85    else:
86      print("This example requires scipy's netcdf support which does not appear to be installed.")
87    

Legend:
Removed from v.4309  
changed lines
  Added in v.4310

  ViewVC Help
Powered by ViewVC 1.1.26