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

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

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

revision 5287 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 magnetic 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    try:
33  try:      import esys.ripley
34     WORKDIR=os.environ['DOWNUNDER_WORKDIR']      HAVE_RIPLEY = True
35    except ImportError:
36        HAVE_RIPLEY = False
37        
38    try:
39        WORKDIR=os.environ['DOWNUNDER_WORKDIR']
40  except KeyError:  except KeyError:
41     WORKDIR='.'      WORKDIR='.'
42    
43  features=[SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
44       x=8*U.km, y=3*U.km, depth=2.5*U.km, v_inner=2., v_outer=1e-6),\  class Test_MagneticInversion2D(unittest.TestCase):
45            SmoothAnomaly(lx=25*U.km, ly=20*U.km, lz=20*U.km,      def test_2D_inversion(self):
46       x=30*U.km, y=1*U.km, depth=18*U.km, v_inner=-1., v_outer=1e-6),\          logging.getLogger('inv.MinimizerLBFGS').setLevel(logging.CRITICAL)
47            SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \          logging.getLogger('inv.MagneticInversion').setLevel(logging.CRITICAL)
48       x=68*U.km, y=3*U.km, depth=5*U.km, v_inner=20., v_outer=1e-6)]          # interesting parameters:
49            depth_offset = 0. * U.km
50  B_b = [2201.*U.Nano*U.Tesla, 31232.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]          n_humps_h = 1
51            n_humps_v = 1
52  source=SyntheticFeatureData(DataSource.MAGNETIC, DIM=2, number_of_elements=30, length=100*U.km, features=features, B_b=B_b)          mu = 1.
53            n_cells_in_data = 200
54  domainbuilder=DomainBuilder(dim=2)          full_knowledge = False
55  domainbuilder.addSource(source)          B_b = [2201.*U.Nano*U.Tesla, 31232.*U.Nano*U.Tesla, -41405.*U.Nano*U.Tesla]
56  domainbuilder.setElementPadding(10)          #
57  domainbuilder.setVerticalExtents(depth=30*U.km, air_layer=10*U.km, num_cells=16)          DIM = 2
58  domainbuilder.setBackgroundMagneticFluxDensity(B_b)          n_cells_in_data = max(n_humps_h*7, n_cells_in_data)
59            l_data = 100. * U.km
60  inv=MagneticInversion()          l_pad = 40. * U.km
61  inv.setSolverTolerance(1e-4)          THICKNESS = 20. * U.km
62  inv.setSolverMaxIterations(10)          l_air = 20. * U.km
63  inv.setup(domainbuilder)          n_cells_v = max(
64  k_new = inv.run()                  int((2*l_air+THICKNESS+depth_offset)/l_data*n_cells_in_data + 0.5), 25)
65  B, chi = inv.getCostFunction().getForwardModel().getSurvey(0)  
66  saveSilo(os.path.join(WORKDIR, 'maginv'), sus=k_new, sus_ref=source.getReferenceProperty(), B=B, chi=chi)  
67            source=SyntheticData(DataSource.MAGNETIC, n_length=n_humps_h, n_depth=n_humps_v,
68                    depth=THICKNESS+depth_offset, depth_offset=depth_offset,
69                    DIM=DIM, number_of_elements=n_cells_in_data, length=l_data, B_b=B_b,
70                    data_offset=0, full_knowledge=full_knowledge)
71    
72            domainbuilder=DomainBuilder(dim=DIM)
73            domainbuilder.addSource(source)
74            domainbuilder.setVerticalExtents(depth=l_air+THICKNESS+depth_offset,
75                                             air_layer=l_air, num_cells=n_cells_v)
76            domainbuilder.setBackgroundMagneticFluxDensity(B_b)
77            domainbuilder.setPadding(l_pad)
78            domainbuilder.fixSusceptibilityBelow(depth=THICKNESS+depth_offset)
79    
80            inv=MagneticInversion()
81            inv.setSolverTolerance(1e-4)
82            inv.setSolverMaxIterations(50)
83            inv.setup(domainbuilder)
84            inv.getCostFunction().setTradeOffFactorsModels(mu)
85    
86            k_new=inv.run()
87            k_ref=source.getReferenceProperty()
88            B, chi = inv.getCostFunction().getForwardModel().getSurvey(0)
89            #saveSilo(os.path.join(WORKDIR, 'results_magnetic_2d'), k=k_new, k_ref=k_ref, B=B, chi=chi)
90    
91    if __name__ == '__main__':
92        run_tests(__name__, exit_on_failure=True)

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26