/[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 4981 - (hide annotations)
Mon Jun 2 00:26:35 2014 UTC (4 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 14205 byte(s)
MT tests added
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 jfenwick 4938 import esys.escriptcore.utestselect as unittest
26 gross 4688 import numpy as np
27     import os
28     import sys
29 gross 4979 import cmath
30 gross 4688 from esys.downunder import *
31     from esys.escript import unitsSI as U
32     from esys.escript import *
33     from esys.weipa import saveSilo
34     from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
35 jfenwick 4713 from esys.escript import getEscriptParamInt
36 gross 4688
37     mpisize = getMPISizeWorld()
38     # this is mainly to avoid warning messages
39     logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
40    
41     try:
42     TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
43     except KeyError:
44     TEST_DATA_ROOT='ref_data'
45    
46     try:
47     WORKDIR=os.environ['DOWNUNDER_WORKDIR']
48     except KeyError:
49     WORKDIR='.'
50 jfenwick 4713
51 gross 4688
52 jfenwick 4713 have_direct=getEscriptParamInt("PASO_DIRECT")
53 gross 4688
54    
55 jfenwick 4713 @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
56 gross 4688 class TestAcousticInversion(unittest.TestCase):
57     def test_API(self):
58     from esys.ripley import Rectangle
59     domain=Rectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
60     omega=2.
61    
62    
63     data=Data([1,2], FunctionOnBoundary(domain))
64     F=Data([2,3], Function(domain))
65     w=1.
66     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, data, 1.) # F is a scalar
67     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, [1,2], F) # data is not Data
68     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Data([1,2], Function(domain)), F) # data is not on boundary
69     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, w, Scalar(1, Function(domain)), F) # data is not of shape (2,)
70     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, [1,2], data, F) # w is not a scalar
71     self.assertRaises(ValueError, AcousticWaveForm, domain, omega, Scalar(1, Function(domain)), data, F) # w is not a scalar
72    
73     # now we do a real one
74     acw=AcousticWaveForm(domain, omega, w, data, F)
75     self.assertEqual(acw.getDomain(), domain)
76     pde=acw.setUpPDE()
77     self.assertIsInstance(pde, LinearPDE)
78     self.assertEqual(pde.getNumEquations(), 2)
79     self.assertEqual(pde.getNumSolutions(), 2)
80     self.assertEqual(pde.getDomain(), domain)
81    
82    
83     def test_numeric2DscaleF(self):
84    
85     from esys.ripley import Rectangle
86     domain=Rectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
87     omega=2.
88    
89     # test solution is u = a * z where a is complex
90     a=complex(3.45, 0.56)
91     sigma=complex(1e-3, 0.056)
92    
93    
94     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
95     mydata=data.copy()
96    
97     z=FunctionOnBoundary(domain).getX()[1]
98     w=whereZero(z-1.)
99     # source:
100     F=Data( [1,0],Function(domain))
101     #
102     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
103     # check rescaled data
104     surv=acw.getSurvey()
105     self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
106    
107     mydata_scale=sqrt( integrate(w*length(mydata)**2) )
108     self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
109     self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
110    
111     # this should be zero:
112     sigma_comps=[sigma.real, sigma.imag]
113     args=acw.getArguments(sigma_comps)
114     d=acw.getDefect(sigma_comps, *args)
115     self.assertTrue(isinstance(d, float))
116     self.assertTrue(abs(d) < 1e-10)
117    
118     dg=acw.getGradient(sigma_comps, *args)
119     self.assertTrue(isinstance(dg, Data))
120     self.assertTrue(dg.getShape()==(2,))
121     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
122     self.assertTrue(Lsup(dg) < 1e-10)
123    
124     # this shuld be zero' too
125     sigma_comps=[2*sigma.real, sigma.imag/2.]
126     args=acw.getArguments(sigma_comps)
127     d=acw.getDefect(sigma_comps, *args)
128     self.assertTrue(isinstance(d, float))
129     self.assertTrue(abs(d)< 1e-10)
130    
131     dg=acw.getGradient(sigma_comps, *args)
132     self.assertTrue(isinstance(dg, Data))
133     self.assertTrue(dg.getShape()==(2,))
134     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
135     self.assertTrue(Lsup(dg) < 1e-10)
136    
137     # this shouldn't be zero:
138     sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
139     args=acw.getArguments(sigma0)
140     d0=acw.getDefect(sigma0, *args)
141     self.assertTrue(isinstance(d0, float))
142     self.assertTrue(d0 >= 0)
143     self.assertTrue(d0 > 1e-10)
144    
145     dg0=acw.getGradient(sigma0, *args)
146     self.assertTrue(isinstance(dg0, Data))
147     self.assertTrue(dg0.getShape()==(2,))
148     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
149     self.assertTrue(Lsup(dg0) > 1e-10)
150    
151     # test the gradient numerrically:
152     h=0.002
153     X=Function(domain).getX()
154     # .. increment:
155     p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
156    
157    
158     sigma1=sigma0+p*[1,0]
159     args=acw.getArguments(sigma1)
160     d1=acw.getDefect(sigma1, *args)
161     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
162    
163     sigma2=sigma0+p*[0,1]
164     args=acw.getArguments(sigma2)
165     d2=acw.getDefect(sigma2, *args)
166     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
167    
168     def test_numeric2DnoscaleF(self):
169    
170     from esys.ripley import Rectangle
171     domain=Rectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
172     omega=1.5
173    
174     # test solution is u = a * z where a is complex
175     a=complex(3.45, 0.56)
176     sigma=complex(1e-3, 0.056)
177    
178    
179     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
180     z=FunctionOnBoundary(domain).getX()[1]
181     w=whereZero(z-1.)
182     # F = - a*omega* sigma
183     F=Data( [-(a*omega**2*sigma).real, -(a*omega**2*sigma).imag ],Function(domain))
184    
185     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=False)
186     # this should be zero:
187     sigma_comps=[sigma.real, sigma.imag]
188     args=acw.getArguments(sigma_comps)
189     d=acw.getDefect(sigma_comps, *args)
190     self.assertTrue(isinstance(d, float))
191 jfenwick 4719 self.assertTrue(Lsup(d) < 1e-10)
192     #self.assertTrue(d >= 0)
193     #self.assertTrue(d < 1e-10)
194 gross 4688
195     dg=acw.getGradient(sigma_comps, *args)
196    
197     self.assertTrue(isinstance(dg, Data))
198     self.assertTrue(dg.getShape()==(2,))
199     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
200     self.assertTrue(Lsup(dg) < 5e-10)
201     # this shouldn't be zero:
202     sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
203     args=acw.getArguments(sigma0)
204     d0=acw.getDefect(sigma0, *args)
205     self.assertTrue(isinstance(d0, float))
206     self.assertTrue(d0 >= 0)
207     self.assertTrue(d0 > 1e-10)
208    
209     dg0=acw.getGradient(sigma0, *args)
210     self.assertTrue(isinstance(dg0, Data))
211     self.assertTrue(dg0.getShape()==(2,))
212     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
213     self.assertTrue(Lsup(dg0) > 1e-10)
214     # test the gradient numerrically:
215     h=0.001
216     X=Function(domain).getX()
217     p=h*sin(length(X)*np.pi)*Lsup(length(sigma0))
218    
219     sigma1=sigma0+p*[1,0]
220     args=acw.getArguments(sigma1)
221     d1=acw.getDefect(sigma1, *args)
222    
223     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
224    
225     sigma2=sigma0+p*[0,1]
226     args=acw.getArguments(sigma2)
227     d2=acw.getDefect(sigma2, *args)
228     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
229    
230 gross 4979 class TestMT2DModelTEMode(unittest.TestCase):
231     def test_API(self):
232     from esys.ripley import Rectangle
233     domain=Rectangle(20,20)
234     omega=2.
235     x=[ [0.2,0.5], [0.3,0.5] ]
236     Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
237     eta=1.
238     w0=1.
239     E_x0=1.
240     # now we do a real one
241     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, E_x0=E_x0)
242     self.assertEqual(acw.getDomain(), domain)
243     pde=acw.setUpPDE()
244     self.assertIsInstance(pde, LinearPDE)
245     self.assertEqual(pde.getNumEquations(), 2)
246     self.assertEqual(pde.getNumSolutions(), 2)
247     self.assertEqual(pde.getDomain(), domain)
248    
249     # other things that should work
250     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
251    
252     # these shouldn't work
253     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
254     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, Z_XY, eta=[1.], w0=[2.,3.], E_x0=complex(4.5,6) )
255     self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, [(6.7,5)], Z_XY, eta=[1.,1.], w0=[2.,3.], E_x0=complex(4.5,6) )
256    
257     def test_PDE(self):
258    
259     omega=2.
260     mu0=0.123
261     SIGMA=15.
262     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
263    
264     from esys.ripley import Rectangle
265     domain=Rectangle(200,200)
266    
267    
268     IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)
269 gross 4980 Z_XY=[ IMP, IMP ]
270 gross 4979 x=[ [0.3,0.5], [0.6,0.5] ]
271 gross 4980 eta=0.005
272 gross 4979 z=domain.getX()[1]
273     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
274     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
275     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
276     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
277    
278     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, E_x0=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.] )
279    
280     args=acw.getArguments(SIGMA)
281     Ex=args[0]
282     Exz=args[1]
283     self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
284     self.assertTrue(Lsup(Ex[1]-Ex1_ex) <= 1e-4 * Lsup(Ex1_ex))
285     self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
286     self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
287    
288 gross 4980 argsr=acw.getArguments(0.)
289     ref=acw.getDefect(0., *argsr)
290    
291     # this should be almost zero:
292 gross 4979 args=acw.getArguments(SIGMA)
293     d=acw.getDefect(SIGMA, *args)
294    
295     self.assertTrue( d > 0.)
296     self.assertTrue( ref > 0.)
297     self.assertTrue( d <= 1e-4 * ref ) # d should be zero (some sort of)
298 gross 4980
299     # and this should be zero
300     d0=acw.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
301     self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
302 gross 4979
303 gross 4980
304     # and this too
305     dg=acw.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
306     self.assertTrue(isinstance(dg, Data))
307     self.assertTrue(dg.getShape()==())
308     self.assertTrue(Lsup(dg) < 1e-10)
309 gross 4979
310 gross 4981 def test_Differential(self):
311 gross 4980
312 gross 4981 INC=0.01
313 gross 4980
314     omega=2.
315     mu0=0.123
316     SIGMA=15.
317     k=cmath.sqrt(1j*omega*mu0*SIGMA) # Ex=exp(k*z)
318    
319     from esys.ripley import Rectangle
320     domain=Rectangle(200,200)
321    
322    
323     IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)
324     Z_XY=[ IMP, IMP ]
325     x=[ [0.3,0.5], [0.6,0.5] ]
326     eta=0.005
327     z=domain.getX()[1]
328     Ex0_ex=exp(-k.real*z)*cos(-k.imag*z)
329     Ex0_ex_z=(-k.real*cos(-k.imag*z)+k.imag*sin(-k.imag*z)) * exp(-k.real*z)
330     Ex1_ex=exp(-k.real*z)*sin(-k.imag*z)
331     Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)
332    
333     acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, E_x0=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.] )
334 gross 4981
335     # this is the base line:
336     SIGMA0=10.
337 gross 4980 args0=acw.getArguments(SIGMA0)
338     d0=acw.getDefect(SIGMA0, *args0)
339    
340 gross 4981 dg0=acw.getGradient(SIGMA0, *args0)
341     self.assertTrue(isinstance(dg0, Data))
342     self.assertTrue(dg0.getShape()==())
343 gross 4980
344 gross 4981
345 gross 4980 X=Function(domain).getX()
346 gross 4981
347     # test 1
348     p=INC
349     SIGMA1=SIGMA0+p
350 gross 4980 args1=acw.getArguments(SIGMA1)
351     d1=acw.getDefect(SIGMA1, *args1)
352 gross 4981 self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
353    
354     # test 2
355     p=exp(-length(X-(0.2,0.2))**2/10)*INC
356     SIGMA1=SIGMA0+p
357     args1=acw.getArguments(SIGMA1)
358     d1=acw.getDefect(SIGMA1, *args1)
359     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
360 gross 4980
361 gross 4981 # test 3
362     p=sin(length(X)*3*3.14)*INC
363     SIGMA1=SIGMA0+p
364     args1=acw.getArguments(SIGMA1)
365     d1=acw.getDefect(SIGMA1, *args1)
366     self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2 * abs(d1-d0) )
367 gross 4979
368 gross 4688
369     if __name__ == "__main__":
370     suite = unittest.TestSuite()
371 gross 4981 suite.addTest(unittest.makeSuite(TestAcousticInversion))
372 gross 4979 suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))
373 gross 4688 s=unittest.TextTestRunner(verbosity=2).run(suite)
374     if not s.wasSuccessful(): sys.exit(1)
375    

  ViewVC Help
Powered by ViewVC 1.1.26