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

Diff of /branches/doubleplusgood/doc/examples/inversion/grav_netcdf.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 26  __url__="https://launchpad.net/escript-f Line 26  __url__="https://launchpad.net/escript-f
26  from esys.downunder import *  from esys.downunder import *
27  from esys.weipa import *  from esys.weipa import *
28  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
29    from esys.escript import saveDataCSV
30    
31    
32    
33    
34  # Set parameters  # Set parameters
35  DATASET = 'data/QLDWest_grav.nc'  DATASET = 'data/GravitySmall.nc'
36    DATA_UNITS = 1e-6 * U.m/(U.sec**2)
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 37  n_cells_v = 25 Line 42  n_cells_v = 25
42  MU = 0.1  MU = 0.1
43    
44  # Setup and run the inversion  # Setup and run the inversion
45  source=NetCdfData(NetCdfData.GRAVITY, DATASET)  def work():
46  db=DomainBuilder(dim=3)    source=NetCdfData(NetCdfData.GRAVITY, DATASET, scale_factor=DATA_UNITS)
47  db.addSource(source)    db=DomainBuilder(dim=3)
48  db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)    db.addSource(source)
49  db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)    db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v)
50  db.fixDensityBelow(depth=thickness)    db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y)
51      db.fixDensityBelow(depth=thickness)
52  inv=GravityInversion()  
53  inv.setSolverTolerance(1e-4)    inv=GravityInversion()
54  inv.setSolverMaxIterations(50)    inv.setSolverTolerance(1e-4)
55  inv.setup(db)    inv.setSolverMaxIterations(50)
56  inv.getCostFunction().setTradeOffFactorsModels(MU)    inv.setup(db)
57      inv.getCostFunction().setTradeOffFactorsModels(MU)
58  print("Starting inversion, please stand by...")  
59  rho = inv.run()    density = inv.run()
60  print("density = %s"%density)    print("density = %s"%density)
61    
62  print("Results saved in result0.silo:")    g, w =  db.getGravitySurveys()[0]
63  g, w =  inv.getCostFunction().getForwardModel().getSurvey(0)    saveSilo("result_gravity.silo", density=density, gravity_anomaly=g, gravity_weight=w)
64  saveSilo("result0.silo", density=rho, gravity_anomaly=g[2], gravity_weight=w[2])    print("Results saved in result_gravity.silo")
65    
66  print("Results saved in result0.vtu:")    saveVTK("result_gravity.vtu", density=density, gravity_anomaly=g, gravity_weight=w)
67  saveSilo("result0.vtu", density=rho)    print("Results saved in result_gravity.vtu")
68    
69  print("Results saved in result0.csv:")    saveDataCSV("result_gravity.csv", density=density, x=density.getFunctionSpace().getX())
70  saveDataCSV("result0.csv", density=rho, x=rho.getFunctionSpace().getX())    print("Results saved in result_gravity.csv")
71    
72      print("All done. Have a nice day!")
73    
74    if 'NetCdfData' in dir():
75      work()
76    else:
77      print("This example requires scipy's netcdf support which does not appear to be installed.")
78    
 print("All done. Have a nice day.!")  

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

  ViewVC Help
Powered by ViewVC 1.1.26