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