| 1 |
|
from __future__ import print_function |
| 2 |
############################################################################## |
############################################################################## |
| 3 |
# |
# |
| 4 |
# Copyright (c) 2003-2013 by University of Queensland |
# Copyright (c) 2003-2015 by The University of Queensland |
| 5 |
# http://www.uq.edu.au |
# http://www.uq.edu.au |
| 6 |
# |
# |
| 7 |
# Primary Business: Queensland, Australia |
# Primary Business: Queensland, Australia |
| 9 |
# http://www.opensource.org/licenses/osl-3.0.php |
# http://www.opensource.org/licenses/osl-3.0.php |
| 10 |
# |
# |
| 11 |
# Development until 2012 by Earth Systems Science Computational Center (ESSCC) |
# Development until 2012 by Earth Systems Science Computational Center (ESSCC) |
| 12 |
# Development since 2012 by School of Earth Sciences |
# Development 2012-2013 by School of Earth Sciences |
| 13 |
|
# Development from 2014 by Centre for Geoscience Computing (GeoComp) |
| 14 |
# |
# |
| 15 |
############################################################################## |
############################################################################## |
| 16 |
|
|
| 17 |
__copyright__="""Copyright (c) 2003-2013 by University of Queensland |
"""2D gravity inversion example using synthetic data""" |
| 18 |
|
|
| 19 |
|
__copyright__="""Copyright (c) 2003-2015 by The University of Queensland |
| 20 |
http://www.uq.edu.au |
http://www.uq.edu.au |
| 21 |
Primary Business: Queensland, Australia""" |
Primary Business: Queensland, Australia""" |
| 22 |
__license__="""Licensed under the Open Software License version 3.0 |
__license__="""Licensed under the Open Software License version 3.0 |
| 23 |
http://www.opensource.org/licenses/osl-3.0.php""" |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 24 |
__url__="https://launchpad.net/escript-finley" |
__url__="https://launchpad.net/escript-finley" |
| 25 |
|
|
|
import logging |
|
| 26 |
import os |
import os |
| 27 |
|
import esys.escriptcore.utestselect as unittest |
| 28 |
|
from esys.escriptcore.testing import * |
| 29 |
from esys.downunder import * |
from esys.downunder import * |
| 30 |
from esys.escript import unitsSI as U |
from esys.escript import unitsSI as U |
| 31 |
from esys.weipa import saveSilo |
from esys.weipa import saveSilo |
| 32 |
|
import logging |
| 33 |
|
|
| 34 |
|
try: |
| 35 |
try: |
import esys.ripley |
| 36 |
WORKDIR=os.environ['DOWNUNDER_WORKDIR'] |
HAVE_RIPLEY = True |
| 37 |
|
except ImportError: |
| 38 |
|
HAVE_RIPLEY = False |
| 39 |
|
|
| 40 |
|
try: |
| 41 |
|
WORKDIR=os.environ['DOWNUNDER_WORKDIR'] |
| 42 |
except KeyError: |
except KeyError: |
| 43 |
WORKDIR='.' |
WORKDIR='.' |
| 44 |
|
|
| 45 |
features=[SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \ |
@unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available") |
| 46 |
x=22*U.km, y=3*U.km, depth=10*U.km, v_inner=200., v_outer=1e-6),\ |
class Test_GravityInversion2D(unittest.TestCase): |
| 47 |
SmoothAnomaly(lx=25*U.km, ly=20*U.km, lz=20*U.km, |
def test_inversion(self): |
| 48 |
x=40*U.km, y=1*U.km, depth=22*U.km, v_inner=-500., v_outer=1e-6),\ |
#be quiet about it |
| 49 |
SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \ |
logging.getLogger('inv.MinimizerLBFGS').setLevel(logging.CRITICAL) |
| 50 |
x=68*U.km, y=3*U.km, depth=13*U.km, v_inner=200., v_outer=1e-6)] |
logging.getLogger('inv.GravityInversion').setLevel(logging.CRITICAL) |
| 51 |
|
|
| 52 |
logger=logging.getLogger('inv') |
# interesting parameters: |
| 53 |
logger.setLevel(logging.DEBUG) |
n_humps_h = 3 |
| 54 |
handler=logging.StreamHandler() |
n_humps_v = 1 |
| 55 |
handler.setLevel(logging.DEBUG) |
mu = 100 |
| 56 |
logger.addHandler(handler) |
n_cells_in_data = 100 |
| 57 |
source=SyntheticFeatureData(DataSource.GRAVITY, DIM=2, number_of_elements=220, length=100*U.km, features=features) |
# ignore: |
| 58 |
domainbuilder=DomainBuilder(dim=2) |
full_knowledge = False |
| 59 |
domainbuilder.addSource(source) |
depth_offset = 0. * U.km |
| 60 |
domainbuilder.setElementPadding(20) |
# |
| 61 |
domainbuilder.setVerticalExtents(depth=50*U.km, air_layer=20*U.km, num_cells=25) |
DIM = 2 |
| 62 |
|
n_cells_in_data = max(n_humps_h*7, n_cells_in_data) |
| 63 |
inv=GravityInversion() |
l_data = 100. * U.km |
| 64 |
inv.setup(domainbuilder) |
l_pad = 40. * U.km |
| 65 |
inv.setSolverTolerance(1e-4) |
THICKNESS = 20. * U.km |
| 66 |
inv.setSolverMaxIterations(40) |
l_air = 20. * U.km |
| 67 |
inv.getCostFunction().setTradeOffFactorsModels(10) |
n_cells_v = max(int((2*l_air+THICKNESS+depth_offset)/ \ |
| 68 |
|
l_data*n_cells_in_data + 0.5), 25) |
| 69 |
rho_new=inv.run() |
|
| 70 |
print "rho_new = ",rho_new |
|
| 71 |
print "rho =", source.getReferenceProperty() |
source=SyntheticData(DataSource.GRAVITY, n_length=n_humps_h, |
| 72 |
g, chi = inv.getCostFunction().getForwardModel().getSurvey(0) |
n_depth=n_humps_v, |
| 73 |
saveSilo(os.path.join(WORKDIR, 'gravinv'), density=rho_new, density_ref=source.getReferenceProperty(), g=g, chi=chi) |
depth=THICKNESS+depth_offset, depth_offset=depth_offset, |
| 74 |
|
DIM=DIM, number_of_elements=n_cells_in_data, length=l_data, |
| 75 |
|
data_offset=0, full_knowledge=full_knowledge) |
| 76 |
|
|
| 77 |
|
domainbuilder=DomainBuilder(dim=DIM) |
| 78 |
|
domainbuilder.addSource(source) |
| 79 |
|
domainbuilder.setVerticalExtents(depth=l_air+THICKNESS+depth_offset, |
| 80 |
|
air_layer=l_air, num_cells=n_cells_v) |
| 81 |
|
domainbuilder.setPadding(l_pad) |
| 82 |
|
domainbuilder.fixDensityBelow(depth=THICKNESS+depth_offset) |
| 83 |
|
|
| 84 |
|
inv=GravityInversion() |
| 85 |
|
inv.setSolverTolerance(1e-4) |
| 86 |
|
inv.setSolverMaxIterations(50) |
| 87 |
|
inv.setup(domainbuilder) |
| 88 |
|
inv.getCostFunction().setTradeOffFactorsModels(mu) |
| 89 |
|
|
| 90 |
|
rho_new = inv.run() |
| 91 |
|
rho_ref = source.getReferenceProperty() |
| 92 |
|
g, chi = inv.getCostFunction().getForwardModel().getSurvey(0) |
| 93 |
|
|
| 94 |
|
if __name__ == '__main__': |
| 95 |
|
run_tests(__name__, exit_on_failure=True) |
Parent Directory
| 
Revision Log
| 
Patch