/[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 3947 by caltinay, Wed Aug 22 23:19:10 2012 UTC revision 5593 by jfenwick, Fri Apr 24 01:36:26 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 The University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
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)
12    # Development 2012-2013 by School of Earth Sciences
13    # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14    #
15    ##############################################################################
16    
17    """2D gravity inversion example using synthetic data"""
18    
19  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
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 *
30  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
31  from esys.escript import inf,sup  from esys.weipa import saveSilo
32  from esys.downunder.datasources import SyntheticDataSource,SmoothAnomaly  import logging
 from esys.downunder.inversions import GravityInversion  
   
33    
34  try:  try:
35     WORKDIR=os.environ['DOWNUNDER_WORKDIR']      import esys.ripley
36        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=50*U.km, ly=20*U.km, lz=40*U.km, \  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
46       x=100*U.km, y=3*U.km, depth=25*U.km, rho_inner=200., rho_outer=1e-6),\  class Test_GravityInversion2D(unittest.TestCase):
47            SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km,      def test_inversion(self):
48       x=400*U.km, y=1*U.km, depth=40*U.km, rho_inner=-200, rho_outer=1e-6)]          #be quiet about it
49            logging.getLogger('inv.MinimizerLBFGS').setLevel(logging.CRITICAL)
50  logger=logging.getLogger('inv')          logging.getLogger('inv.GravityInversion').setLevel(logging.CRITICAL)
51  logger.setLevel(logging.FATAL)  
52  handler=logging.StreamHandler()          # interesting parameters:
53  handler.setLevel(logging.FATAL)          n_humps_h = 3
54  logger.addHandler(handler)          n_humps_v = 1
55  source=SyntheticDataSource(DIM=2, NE=20, l=500*U.km, h=60*U.km, features=features)          mu = 100
56  source.setPadding(5, 0.1)          n_cells_in_data = 100
57  inv=GravityInversion()          # ignore:
58  inv.setDataSource(source)          full_knowledge = False
59  inv.setOutputDirectory(WORKDIR)          depth_offset = 0. * U.km
60  inv.setSolverTolerance(1e-5)          #
61  inv.setSolverMaxIterations(100)          DIM = 2
62  inv.setSolverOptions(initialHessian=100)          n_cells_in_data = max(n_humps_h*7, n_cells_in_data)
63  x=source.getDomain().getX()          l_data = 100. * U.km
64  l0=sup(x[0])-inf(x[0])          l_pad = 40. * U.km
65  l1=sup(x[1])-inf(x[1])          THICKNESS = 20. * U.km
66  l=max(l0,l1)          l_air = 20. * U.km
67  G=6.6742e-11          n_cells_v = max(int((2*l_air+THICKNESS+depth_offset)/ \
68  mu=0.5*(l**2*G)**2                          l_data*n_cells_in_data + 0.5), 25)
69  inv.setWeights(mu_reg=mu)  
70  rho_new=inv.run()  
71            source=SyntheticData(DataSource.GRAVITY, n_length=n_humps_h,
72                    n_depth=n_humps_v,
73                    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.3947  
changed lines
  Added in v.5593

  ViewVC Help
Powered by ViewVC 1.1.26