| 1 |
gross |
4116 |
|
| 2 |
|
|
############################################################################## |
| 3 |
|
|
# |
| 4 |
|
|
# Copyright (c) 2003-2012 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 |
|
|
__copyright__="""Copyright (c) 2003-2012 by University of Queensland |
| 17 |
|
|
http://www.uq.edu.au |
| 18 |
|
|
Primary Business: Queensland, Australia""" |
| 19 |
|
|
__license__="""Licensed under the Open Software License version 3.0 |
| 20 |
|
|
http://www.opensource.org/licenses/osl-3.0.php""" |
| 21 |
|
|
__url__="https://launchpad.net/escript-finley" |
| 22 |
|
|
|
| 23 |
|
|
import logging |
| 24 |
|
|
import os |
| 25 |
|
|
from esys.downunder import * |
| 26 |
|
|
from esys.escript import unitsSI as U |
| 27 |
|
|
from esys.weipa import saveSilo |
| 28 |
|
|
|
| 29 |
|
|
|
| 30 |
|
|
try: |
| 31 |
|
|
WORKDIR=os.environ['DOWNUNDER_WORKDIR'] |
| 32 |
|
|
except KeyError: |
| 33 |
|
|
WORKDIR='.' |
| 34 |
|
|
logger=logging.getLogger('inv') |
| 35 |
|
|
logger.setLevel(logging.DEBUG) |
| 36 |
|
|
handler=logging.StreamHandler() |
| 37 |
|
|
handler.setLevel(logging.DEBUG) |
| 38 |
|
|
logger.addHandler(handler) |
| 39 |
|
|
|
| 40 |
|
|
# interesting parameter: |
| 41 |
gross |
4120 |
depth_offset=0.*U.km |
| 42 |
|
|
n_humbs_h= 3 |
| 43 |
|
|
n_humbs_v=1 |
| 44 |
|
|
mu=100 |
| 45 |
|
|
n_cells_in_data=100 |
| 46 |
gross |
4116 |
full_knowledge=False |
| 47 |
|
|
# |
| 48 |
|
|
n_cells_in_data=max(n_humbs_h*7,n_cells_in_data) |
| 49 |
|
|
l_data = 100 * U.km |
| 50 |
|
|
l_pad=40*U.km |
| 51 |
|
|
THICKNESS=20.*U.km |
| 52 |
|
|
l_air=20*U.km |
| 53 |
|
|
n_cells_v=max(int((2*l_air+THICKNESS+depth_offset)/l_data*n_cells_in_data + 0.5), 25) |
| 54 |
|
|
|
| 55 |
|
|
|
| 56 |
|
|
source=SyntheticData(DataSource.GRAVITY,n_length=n_humbs_h, n_depth=n_humbs_v, depth=THICKNESS+depth_offset, depth_offset=depth_offset, |
| 57 |
|
|
DIM=2, number_of_elements =n_cells_in_data, length=l_data, |
| 58 |
|
|
data_offset=0,full_knowledge=full_knowledge, spherical=False) |
| 59 |
|
|
|
| 60 |
|
|
|
| 61 |
|
|
domainbuilder=DomainBuilder(dim=2) |
| 62 |
|
|
domainbuilder.addSource(source) |
| 63 |
|
|
domainbuilder.setVerticalExtents(depth=l_air+THICKNESS+depth_offset, air_layer=l_air, num_cells=n_cells_v) |
| 64 |
|
|
domainbuilder.setPadding(l_pad) |
| 65 |
gross |
4120 |
domainbuilder.fixDensityBelow(depth=THICKNESS+depth_offset) |
| 66 |
gross |
4116 |
|
| 67 |
|
|
inv=GravityInversion() |
| 68 |
|
|
inv.setSolverTolerance(1e-4) |
| 69 |
|
|
inv.setSolverMaxIterations(50) |
| 70 |
|
|
inv.setup(domainbuilder) |
| 71 |
gross |
4122 |
inv.getCostFunction().setTradeOffFactorsModels(mu) |
| 72 |
gross |
4116 |
|
| 73 |
|
|
rho_new=inv.run() |
| 74 |
|
|
print "rho_new = ",rho_new |
| 75 |
|
|
print "rho =", source.getReferenceProperty() |
| 76 |
gross |
4122 |
g, chi = inv.getCostFunction().getForwardModels()[0].getSurvey(0) |
| 77 |
gross |
4116 |
saveSilo(os.path.join(WORKDIR, 'gravinv'), density=rho_new, density_ref=source.getReferenceProperty(), g=g, chi=chi) |
| 78 |
|
|
|
Parent Directory
| 
Revision Log