/[escript]/trunk/downunder/test/python/run_gravity.py
ViewVC logotype

Diff of /trunk/downunder/test/python/run_gravity.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4115 by gross, Fri Dec 14 04:48:48 2012 UTC revision 5448 by jfenwick, Fri Feb 6 05:31:37 2015 UTC
# Line 1  Line 1 
1    from __future__ import print_function
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2015 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 9  Line 9 
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-2012 by University of Queensland  """2D gravity inversion example using synthetic data"""
18    
19    __copyright__="""Copyright (c) 2003-2015 by 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.setSolverTolerance(1e-4)          l_pad = 40. * U.km
65  inv.setSolverMaxIterations(10)          THICKNESS = 20. * U.km
66  #inv.setTradeOffFactors(mu_model=30.)          l_air = 20. * U.km
67  inv.setup(domainbuilder)          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.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)

Legend:
Removed from v.4115  
changed lines
  Added in v.5448

  ViewVC Help
Powered by ViewVC 1.1.26