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

Annotation of /trunk/downunder/test/python/run_forward.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4713 - (hide annotations)
Thu Feb 27 02:00:15 2014 UTC (5 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 9009 byte(s)
Added a temporay fix to get some of the unit tests passing again.
This mechanism will be removed before the release so don't get used to it.
1 gross 4688
2     ##############################################################################
3     #
4     # Copyright (c) 2003-2014 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     __copyright__="""Copyright (c) 2003-2014 by University of Queensland
18     http://www.uq.edu.au
19     Primary Business: Queensland, Australia"""
20     __license__="""Licensed under the Open Software License version 3.0
21     http://www.opensource.org/licenses/osl-3.0.php"""
22     __url__="https://launchpad.net/escript-finley"
23    
24     import logging
25     import unittest
26     import numpy as np
27     import os
28     import sys
29     from esys.downunder import *
30     from esys.escript import unitsSI as U
31     from esys.escript import *
32     from esys.weipa import saveSilo
33     from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
34 jfenwick 4713 from esys.escript import getEscriptParamInt
35 gross 4688
36     mpisize = getMPISizeWorld()
37     # this is mainly to avoid warning messages
38     logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
39    
40     try:
41     TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
42     except KeyError:
43     TEST_DATA_ROOT='ref_data'
44    
45     try:
46     WORKDIR=os.environ['DOWNUNDER_WORKDIR']
47     except KeyError:
48     WORKDIR='.'
49 jfenwick 4713
50 gross 4688
51 jfenwick 4713 have_direct=getEscriptParamInt("PASO_DIRECT")
52 gross 4688
53    
54 jfenwick 4713 @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
55 gross 4688 class TestAcousticInversion(unittest.TestCase):
56     def test_API(self):
57     from esys.ripley import Rectangle
58     domain=Rectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
59     omega=2.
60    
61    
62     data=Data([1,2], FunctionOnBoundary(domain))
63     F=Data([2,3], Function(domain))
64     w=1.
65     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, data, 1.) # F is a scalar
66     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, [1,2], F) # data is not Data
67     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Data([1,2], Function(domain)), F) # data is not on boundary
68     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Scalar(1, Function(domain)), F) # data is not of shape (2,)
69     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, [1,2], data, F) # w is not a scalar
70     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, Scalar(1, Function(domain)), data, F) # w is not a scalar
71    
72     # now we do a real one
73     acw=AcousticWaveForm(domain, omega, w, data, F)
74     self.assertEqual(acw.getDomain(), domain)
75     pde=acw.setUpPDE()
76     self.assertIsInstance(pde, LinearPDE)
77     self.assertEqual(pde.getNumEquations(), 2)
78     self.assertEqual(pde.getNumSolutions(), 2)
79     self.assertEqual(pde.getDomain(), domain)
80    
81    
82     def test_numeric2DscaleF(self):
83    
84     from esys.ripley import Rectangle
85     domain=Rectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
86     omega=2.
87    
88     # test solution is u = a * z where a is complex
89     a=complex(3.45, 0.56)
90     sigma=complex(1e-3, 0.056)
91    
92    
93     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
94     mydata=data.copy()
95    
96     z=FunctionOnBoundary(domain).getX()[1]
97     w=whereZero(z-1.)
98     # source:
99     F=Data( [1,0],Function(domain))
100     #
101     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
102     # check rescaled data
103     surv=acw.getSurvey()
104     self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
105    
106     mydata_scale=sqrt( integrate(w*length(mydata)**2) )
107     self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
108     self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
109    
110     # this should be zero:
111     sigma_comps=[sigma.real, sigma.imag]
112     args=acw.getArguments(sigma_comps)
113     d=acw.getDefect(sigma_comps, *args)
114     self.assertTrue(isinstance(d, float))
115     self.assertTrue(abs(d) < 1e-10)
116    
117     dg=acw.getGradient(sigma_comps, *args)
118     self.assertTrue(isinstance(dg, Data))
119     self.assertTrue(dg.getShape()==(2,))
120     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
121     self.assertTrue(Lsup(dg) < 1e-10)
122    
123     # this shuld be zero' too
124     sigma_comps=[2*sigma.real, sigma.imag/2.]
125     args=acw.getArguments(sigma_comps)
126     d=acw.getDefect(sigma_comps, *args)
127     self.assertTrue(isinstance(d, float))
128     self.assertTrue(abs(d)< 1e-10)
129    
130     dg=acw.getGradient(sigma_comps, *args)
131     self.assertTrue(isinstance(dg, Data))
132     self.assertTrue(dg.getShape()==(2,))
133     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
134     self.assertTrue(Lsup(dg) < 1e-10)
135    
136     # this shouldn't be zero:
137     sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
138     args=acw.getArguments(sigma0)
139     d0=acw.getDefect(sigma0, *args)
140     self.assertTrue(isinstance(d0, float))
141     self.assertTrue(d0 >= 0)
142     self.assertTrue(d0 > 1e-10)
143    
144     dg0=acw.getGradient(sigma0, *args)
145     self.assertTrue(isinstance(dg0, Data))
146     self.assertTrue(dg0.getShape()==(2,))
147     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
148     self.assertTrue(Lsup(dg0) > 1e-10)
149    
150     # test the gradient numerrically:
151     h=0.002
152     X=Function(domain).getX()
153     # .. increment:
154     p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
155    
156    
157     sigma1=sigma0+p*[1,0]
158     args=acw.getArguments(sigma1)
159     d1=acw.getDefect(sigma1, *args)
160     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
161    
162     sigma2=sigma0+p*[0,1]
163     args=acw.getArguments(sigma2)
164     d2=acw.getDefect(sigma2, *args)
165     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
166    
167     def test_numeric2DnoscaleF(self):
168    
169     from esys.ripley import Rectangle
170     domain=Rectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
171     omega=1.5
172    
173     # test solution is u = a * z where a is complex
174     a=complex(3.45, 0.56)
175     sigma=complex(1e-3, 0.056)
176    
177    
178     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
179     z=FunctionOnBoundary(domain).getX()[1]
180     w=whereZero(z-1.)
181     # F = - a*omega* sigma
182     F=Data( [-(a*omega**2*sigma).real, -(a*omega**2*sigma).imag ],Function(domain))
183    
184     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=False)
185     # this should be zero:
186     sigma_comps=[sigma.real, sigma.imag]
187     args=acw.getArguments(sigma_comps)
188     d=acw.getDefect(sigma_comps, *args)
189     self.assertTrue(isinstance(d, float))
190     self.assertTrue(d >= 0)
191     self.assertTrue(d < 1e-10)
192    
193     dg=acw.getGradient(sigma_comps, *args)
194    
195     self.assertTrue(isinstance(dg, Data))
196     self.assertTrue(dg.getShape()==(2,))
197     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
198     self.assertTrue(Lsup(dg) < 5e-10)
199     # this shouldn't be zero:
200     sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
201     args=acw.getArguments(sigma0)
202     d0=acw.getDefect(sigma0, *args)
203     self.assertTrue(isinstance(d0, float))
204     self.assertTrue(d0 >= 0)
205     self.assertTrue(d0 > 1e-10)
206    
207     dg0=acw.getGradient(sigma0, *args)
208     self.assertTrue(isinstance(dg0, Data))
209     self.assertTrue(dg0.getShape()==(2,))
210     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
211     self.assertTrue(Lsup(dg0) > 1e-10)
212     # test the gradient numerrically:
213     h=0.001
214     X=Function(domain).getX()
215     p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
216    
217     sigma1=sigma0+p*[1,0]
218     args=acw.getArguments(sigma1)
219     d1=acw.getDefect(sigma1, *args)
220    
221     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
222    
223     sigma2=sigma0+p*[0,1]
224     args=acw.getArguments(sigma2)
225     d2=acw.getDefect(sigma2, *args)
226     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
227    
228    
229     if __name__ == "__main__":
230     suite = unittest.TestSuite()
231     suite.addTest(unittest.makeSuite(TestAcousticInversion))
232     s=unittest.TextTestRunner(verbosity=2).run(suite)
233     if not s.wasSuccessful(): sys.exit(1)
234    

  ViewVC Help
Powered by ViewVC 1.1.26