| 1 |
gross |
4532 |
|
| 2 |
|
|
############################################################################## |
| 3 |
|
|
# |
| 4 |
|
|
# Copyright (c) 2009-2013 by University of Queensland |
| 5 |
|
|
# http://www.uq.edu.au |
| 6 |
|
|
# |
| 7 |
|
|
# Primary Business: Queensland, Australia |
| 8 |
|
|
# Licensed under the Open Software License version 3.0 |
| 9 |
|
|
# http://www.opensource.org/licenses/osl-3.0.php |
| 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 |
|
|
http://www.opensource.org/licenses/osl-3.0.php""" |
| 23 |
|
|
__url__="https://launchpad.net/escript-finley" |
| 24 |
|
|
|
| 25 |
|
|
# Import required modules |
| 26 |
|
|
from esys.downunder import * |
| 27 |
|
|
from esys.escript import unitsSI as U |
| 28 |
|
|
from esys.escript import saveDataCSV |
| 29 |
|
|
from esys.weipa import * |
| 30 |
|
|
|
| 31 |
|
|
# Set parameters |
| 32 |
|
|
MAGNETIC_DATASET = 'data/MagneticSmall.nc' |
| 33 |
|
|
MAG_UNITS = U.Nano * U.Tesla |
| 34 |
|
|
GRAVITY_DATASET = 'data/GravitySmall.nc' |
| 35 |
|
|
GRAV_UNITS = 1e-6 * U.m/(U.sec**2) |
| 36 |
|
|
# background magnetic field components (B_East, B_North, B_Vertical) |
| 37 |
|
|
B_b = [2201.*U.Nano*U.Tesla, 31232.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla] |
| 38 |
|
|
PAD_X = 0.2 |
| 39 |
|
|
PAD_Y = 0.2 |
| 40 |
|
|
thickness = 40. * U.km |
| 41 |
|
|
l_air = 6. * U.km |
| 42 |
|
|
n_cells_v = 25 |
| 43 |
gross |
4535 |
mu_gravity = 1. |
| 44 |
gross |
4532 |
mu_magnetic = 0.1 |
| 45 |
|
|
COORDINATES=CartesianReferenceSystem() |
| 46 |
|
|
#COORDINATES=WGS84ReferenceSystem() |
| 47 |
|
|
|
| 48 |
|
|
|
| 49 |
|
|
def work(): |
| 50 |
|
|
# Setup and run the inversion |
| 51 |
|
|
grav_source=NetCdfData(NetCdfData.GRAVITY, GRAVITY_DATASET, scale_factor=GRAV_UNITS, reference_system=COORDINATES) |
| 52 |
|
|
mag_source=NetCdfData(NetCdfData.MAGNETIC, MAGNETIC_DATASET, scale_factor=MAG_UNITS, reference_system=COORDINATES) |
| 53 |
|
|
db=DomainBuilder(dim=3, reference_system=COORDINATES) |
| 54 |
|
|
db.addSource(grav_source) |
| 55 |
|
|
db.addSource(mag_source) |
| 56 |
|
|
db.setVerticalExtents(depth=thickness, air_layer=l_air, num_cells=n_cells_v) |
| 57 |
|
|
db.setFractionalPadding(pad_x=PAD_X, pad_y=PAD_Y) |
| 58 |
|
|
db.setBackgroundMagneticFluxDensity(B_b) |
| 59 |
|
|
db.fixDensityBelow(depth=thickness) |
| 60 |
|
|
db.fixSusceptibilityBelow(depth=thickness) |
| 61 |
|
|
|
| 62 |
|
|
inv=StrongJointGravityMagneticInversion() |
| 63 |
|
|
inv.setSolverTolerance(1e-4) |
| 64 |
|
|
inv.setSolverMaxIterations(50) |
| 65 |
|
|
inv.setup(db) |
| 66 |
|
|
inv.getCostFunction().setTradeOffFactorsModels([mu_gravity, mu_magnetic]) |
| 67 |
|
|
inv.getCostFunction().setTradeOffFactorsRegularization(mu = 1.) |
| 68 |
|
|
|
| 69 |
|
|
density, susceptibility = inv.run() |
| 70 |
|
|
print("density = %s"%density) |
| 71 |
|
|
print("susceptibility = %s"%susceptibility) |
| 72 |
|
|
|
| 73 |
|
|
g, wg = db.getGravitySurveys()[0] |
| 74 |
|
|
B, wB = db.getMagneticSurveys()[0] |
| 75 |
caltinay |
4561 |
if saveSilo("result_gravmag_strong.silo", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B, magnetic_weight=wB): |
| 76 |
|
|
print("Results saved in result_gravmag_strong.silo") |
| 77 |
gross |
4532 |
else: |
| 78 |
caltinay |
4561 |
print("Failed to save silo file. Possibly no Silo support.") |
| 79 |
gross |
4532 |
|
| 80 |
caltinay |
4561 |
saveVTK("result_gravmag_strong.vtu", density=density, gravity_anomaly=g, gravity_weight=wg, susceptibility=susceptibility, magnetic_anomaly=B, magnetic_weight=wB) |
| 81 |
|
|
print("Results saved in result_gravmag_strong.vtu") |
| 82 |
gross |
4532 |
|
| 83 |
|
|
print("All done. Have a nice day!") |
| 84 |
|
|
|
| 85 |
|
|
if 'NetCdfData' in dir(): |
| 86 |
|
|
work() |
| 87 |
|
|
else: |
| 88 |
|
|
print("This example requires scipy's netcdf support which does not appear to be installed.") |
| 89 |
|
|
|
Parent Directory
| 
Revision Log