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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5448 - (show annotations)
Fri Feb 6 05:31:37 2015 UTC (4 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 3322 byte(s)
Updating all the dates
1 from __future__ import print_function
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 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 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-2015 by University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 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
31 from esys.weipa import saveSilo
32 import logging
33
34 try:
35 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:
43 WORKDIR='.'
44
45 @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
46 class Test_GravityInversion2D(unittest.TestCase):
47 def test_inversion(self):
48 #be quiet about it
49 logging.getLogger('inv.MinimizerLBFGS').setLevel(logging.CRITICAL)
50 logging.getLogger('inv.GravityInversion').setLevel(logging.CRITICAL)
51
52 # interesting parameters:
53 n_humps_h = 3
54 n_humps_v = 1
55 mu = 100
56 n_cells_in_data = 100
57 # ignore:
58 full_knowledge = False
59 depth_offset = 0. * U.km
60 #
61 DIM = 2
62 n_cells_in_data = max(n_humps_h*7, n_cells_in_data)
63 l_data = 100. * U.km
64 l_pad = 40. * U.km
65 THICKNESS = 20. * U.km
66 l_air = 20. * U.km
67 n_cells_v = max(int((2*l_air+THICKNESS+depth_offset)/ \
68 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)

  ViewVC Help
Powered by ViewVC 1.1.26