/[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 4980 - (hide annotations)
Sun Jun 1 23:53:11 2014 UTC (4 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 16955 byte(s)
more fixes on MT
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 4980 def ttest_Differential(self):
311    
312     INC=0.001
313    
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    
335     # this should be small
336     args0=acw.getArguments(SIGMA0)
337     d0=acw.getDefect(SIGMA0, *args0)
338     print "d0 =",d0
339    
340     dg=acw.getGradient(SIGMA0, *args0)
341     self.assertTrue(isinstance(dg, Data))
342     self.assertTrue(dg.getShape()==())
343     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
344     print dg
345     self.assertTrue(Lsup(dg) < 1e-10)
346    
347     X=Function(domain).getX()
348     inc=exp(-length(X-(0.2,0.2))**2/0.03)*INC
349    
350     SIGMA1=SIGMA0+inc
351     args1=acw.getArguments(SIGMA1)
352     d1=acw.getDefect(SIGMA1, *args1)
353    
354    
355     print "d0, d1 = ",d0, d1
356    
357    
358    
359    
360 gross 4979 1/0
361     # test solution is u = a * z where a is complex
362     a=complex(3.45, 0.56)
363     sigma=complex(1e-3, 0.056)
364    
365    
366     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
367     mydata=data.copy()
368    
369     z=FunctionOnBoundary(domain).getX()[1]
370     w=whereZero(z-1.)
371     # source:
372     F=Data( [1,0],Function(domain))
373     #
374     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
375     # check rescaled data
376     surv=acw.getSurvey()
377     self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
378    
379     mydata_scale=sqrt( integrate(w*length(mydata)**2) )
380     self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
381     self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
382    
383     # this should be zero:
384     sigma_comps=[sigma.real, sigma.imag]
385     args=acw.getArguments(sigma_comps)
386     d=acw.getDefect(sigma_comps, *args)
387     self.assertTrue(isinstance(d, float))
388     self.assertTrue(abs(d) < 1e-10)
389    
390     dg=acw.getGradient(sigma_comps, *args)
391     self.assertTrue(isinstance(dg, Data))
392     self.assertTrue(dg.getShape()==(2,))
393     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
394     self.assertTrue(Lsup(dg) < 1e-10)
395    
396     # this shuld be zero' too
397     sigma_comps=[2*sigma.real, sigma.imag/2.]
398     args=acw.getArguments(sigma_comps)
399     d=acw.getDefect(sigma_comps, *args)
400     self.assertTrue(isinstance(d, float))
401     self.assertTrue(abs(d)< 1e-10)
402    
403     dg=acw.getGradient(sigma_comps, *args)
404     self.assertTrue(isinstance(dg, Data))
405     self.assertTrue(dg.getShape()==(2,))
406     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
407     self.assertTrue(Lsup(dg) < 1e-10)
408    
409     # this shouldn't be zero:
410     sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
411     args=acw.getArguments(sigma0)
412     d0=acw.getDefect(sigma0, *args)
413     self.assertTrue(isinstance(d0, float))
414     self.assertTrue(d0 >= 0)
415     self.assertTrue(d0 > 1e-10)
416    
417     dg0=acw.getGradient(sigma0, *args)
418     self.assertTrue(isinstance(dg0, Data))
419     self.assertTrue(dg0.getShape()==(2,))
420     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
421     self.assertTrue(Lsup(dg0) > 1e-10)
422    
423     # test the gradient numerrically:
424     h=0.002
425     X=Function(domain).getX()
426     # .. increment:
427     p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
428    
429    
430     sigma1=sigma0+p*[1,0]
431     args=acw.getArguments(sigma1)
432     d1=acw.getDefect(sigma1, *args)
433     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
434    
435     sigma2=sigma0+p*[0,1]
436     args=acw.getArguments(sigma2)
437     d2=acw.getDefect(sigma2, *args)
438     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
439    
440 gross 4688
441     if __name__ == "__main__":
442     suite = unittest.TestSuite()
443 gross 4979 #suite.addTest(unittest.makeSuite(TestAcousticInversion))
444     suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))
445 gross 4688 s=unittest.TextTestRunner(verbosity=2).run(suite)
446     if not s.wasSuccessful(): sys.exit(1)
447    

  ViewVC Help
Powered by ViewVC 1.1.26