# Contents of /trunk/doc/examples/inversion/gravmag_netcdf.py

Revision 4175 - (show annotations)
Wed Jan 30 05:01:50 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 2708 byte(s)
```joint netcdf example and some minor fixes to other example/test case.

```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2009-2013 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 """3D gravity/magnetic joint inversion example using netCDF data""" 17 18 __copyright__="""Copyright (c) 2009-2013 by University of Queensland 19 http://www.uq.edu.au 20 Primary Business: Queensland, Australia""" 21 __license__="""Licensed under the Open Software License version 3.0 22 23 __url__= 24 25 # Import required modules 26 from esys.downunder import * 27 from esys.downunder.minimizers import MinimizerMaxIterReached 28 from esys.escript import unitsSI as U 29 from esys.weipa import saveSilo 30 31 # Set parameters 32 MAGNETIC_DATASET = 'magnetic_anomaly.nc' 33 GRAVITY_DATASET = 'gravity_anomaly.nc' 34 latitude = -28.5 35 PAD_X = 0.2 36 PAD_Y = 0.2 37 thickness = 40. * U.km 38 l_air = 6. * U.km 39 n_cells_v = 25 40 mu_gravity = 10. 41 mu_magnetic = 0.1 42 43 # Setup and run the inversion 44 B_b=simpleGeoMagneticFluxDensity(latitude=latitude) 45 grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET) 46 mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET) 47 db=DomainBuilder(dim=3) 48 db.addSource(grav_source) 49 db.addSource(mag_source) 50 db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v) 51 db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y) 52 db.setBackgroundMagneticFluxDensity(B_b) 53 db.fixDensityBelow(depth=thickness) 54 db.fixSusceptibilityBelow(depth=thickness) 55 56 inv=MagneticInversion() 57 inv.setSolverTolerance(1e-4) 58 inv.setSolverMaxIterations(50) 59 inv.setup(db) 60 inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic]) 61 inv.getCostFunction().setTradeOffFactorsRegularization(mu = [1.,1.], mu_c=1.) 62 63 print("Starting inversion, please stand by...") 64 try: 65 density, susceptibility = inv.run() 66 except MinimizerMaxIterReached as e: 67 print(e) 68 density, susceptibility = inv.p 69 70 print("density = %s"%density) 71 print("susceptibility = %s"%susceptibility) 72 73 # Save results 74 g, wg = inv.getCostFunction().getForwardModel().getSurvey(0) 75 B, wB = inv.getCostFunction().getForwardModel().getSurvey(1) 76 saveSilo("result_gravmag.silo", density=density, gravity_anomaly=g[2], gravity_weight=wg[2], susceptibility=susceptibility, magnetic_anomaly=B[2], magnetic_weight=wB[2]) 77 print("Results saved in result_gravmag.silo. Good bye.") 78

 ViewVC Help Powered by ViewVC 1.1.26