/[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 4657 by jfenwick, Thu Feb 6 06:12:20 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1 
1    from __future__ import print_function
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2014 by University of Queensland  # Copyright (c) 2003-2014 by University of Queensland
# Line 14  Line 14 
14  #  #
15  ##############################################################################  ##############################################################################
16    
17    """2D gravity inversion example using synthetic data"""
18    
19  __copyright__="""Copyright (c) 2003-2014 by University of Queensland  __copyright__="""Copyright (c) 2003-2014 by University of Queensland
20  http://www.uq.edu.au  http://www.uq.edu.au
21  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 22  http://www.opensource.org/licenses/osl-3 Line 24  http://www.opensource.org/licenses/osl-3
24  __url__="https://launchpad.net/escript-finley"  __url__="https://launchpad.net/escript-finley"
25    
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  source=SyntheticFeatureData(DataSource.GRAVITY, DIM=2, number_of_elements=220, length=100*U.km, features=features)          # interesting parameters:
53  domainbuilder=DomainBuilder(dim=2)          n_humps_h = 3
54  domainbuilder.addSource(source)          n_humps_v = 1
55  domainbuilder.setElementPadding(20)          mu = 100
56  domainbuilder.setVerticalExtents(depth=50*U.km, air_layer=20*U.km, num_cells=25)          n_cells_in_data = 100
57            # ignore:
58  inv=GravityInversion()          full_knowledge = False
59  inv.setup(domainbuilder)          depth_offset = 0. * U.km
60  inv.setSolverTolerance(1e-4)          #
61  inv.setSolverMaxIterations(40)          DIM = 2
62  inv.getCostFunction().setTradeOffFactorsModels(10)          n_cells_in_data = max(n_humps_h*7, n_cells_in_data)
63            l_data = 100. * U.km
64  rho_new=inv.run()          l_pad = 40. * U.km
65  print("rho_new = %s"%rho_new)          THICKNESS = 20. * U.km
66  print("rho = %s"%source.getReferenceProperty())          l_air = 20. * U.km
67  g, chi = inv.getCostFunction().getForwardModel().getSurvey(0)          n_cells_v = max(int((2*l_air+THICKNESS+depth_offset)/ \
68  saveSilo(os.path.join(WORKDIR, 'gravinv'), density=rho_new, density_ref=source.getReferenceProperty(), g=g, chi=chi)                          l_data*n_cells_in_data + 0.5), 25)
69    
70    
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.4657  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26