/[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 4979 - (hide annotations)
Sun Jun 1 22:36:16 2014 UTC (4 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 15029 byte(s)
some fixes to the MT forward model
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     Z_XY=[ complex(IMP,0.), complex(IMP,0.) ]
270     x=[ [0.3,0.5], [0.6,0.5] ]
271     eta=0.0005
272     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     args=acw.getArguments(0.)
289     ref=acw.getDefect(0., *args)
290    
291    
292     def ttest_Differential(self):
293    
294    
295     #######
296     args=acw.getArguments(SIGMA)
297     d=acw.getDefect(SIGMA, *args)
298    
299     self.assertTrue( d > 0.)
300     self.assertTrue( ref > 0.)
301     self.assertTrue( d <= 1e-4 * ref ) # d should be zero (some sort of)
302    
303    
304     1/0
305     # test solution is u = a * z where a is complex
306     a=complex(3.45, 0.56)
307     sigma=complex(1e-3, 0.056)
308    
309    
310     data=Data([a.real, a.imag], FunctionOnBoundary(domain))
311     mydata=data.copy()
312    
313     z=FunctionOnBoundary(domain).getX()[1]
314     w=whereZero(z-1.)
315     # source:
316     F=Data( [1,0],Function(domain))
317     #
318     acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)
319     # check rescaled data
320     surv=acw.getSurvey()
321     self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)
322    
323     mydata_scale=sqrt( integrate(w*length(mydata)**2) )
324     self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )
325     self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )
326    
327     # this should be zero:
328     sigma_comps=[sigma.real, sigma.imag]
329     args=acw.getArguments(sigma_comps)
330     d=acw.getDefect(sigma_comps, *args)
331     self.assertTrue(isinstance(d, float))
332     self.assertTrue(abs(d) < 1e-10)
333    
334     dg=acw.getGradient(sigma_comps, *args)
335     self.assertTrue(isinstance(dg, Data))
336     self.assertTrue(dg.getShape()==(2,))
337     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
338     self.assertTrue(Lsup(dg) < 1e-10)
339    
340     # this shuld be zero' too
341     sigma_comps=[2*sigma.real, sigma.imag/2.]
342     args=acw.getArguments(sigma_comps)
343     d=acw.getDefect(sigma_comps, *args)
344     self.assertTrue(isinstance(d, float))
345     self.assertTrue(abs(d)< 1e-10)
346    
347     dg=acw.getGradient(sigma_comps, *args)
348     self.assertTrue(isinstance(dg, Data))
349     self.assertTrue(dg.getShape()==(2,))
350     self.assertTrue(dg.getFunctionSpace()==Solution(domain))
351     self.assertTrue(Lsup(dg) < 1e-10)
352    
353     # this shouldn't be zero:
354     sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
355     args=acw.getArguments(sigma0)
356     d0=acw.getDefect(sigma0, *args)
357     self.assertTrue(isinstance(d0, float))
358     self.assertTrue(d0 >= 0)
359     self.assertTrue(d0 > 1e-10)
360    
361     dg0=acw.getGradient(sigma0, *args)
362     self.assertTrue(isinstance(dg0, Data))
363     self.assertTrue(dg0.getShape()==(2,))
364     self.assertTrue(dg0.getFunctionSpace()==Solution(domain))
365     self.assertTrue(Lsup(dg0) > 1e-10)
366    
367     # test the gradient numerrically:
368     h=0.002
369     X=Function(domain).getX()
370     # .. increment:
371     p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))
372    
373    
374     sigma1=sigma0+p*[1,0]
375     args=acw.getArguments(sigma1)
376     d1=acw.getDefect(sigma1, *args)
377     self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2 * abs(d1-d0) )
378    
379     sigma2=sigma0+p*[0,1]
380     args=acw.getArguments(sigma2)
381     d2=acw.getDefect(sigma2, *args)
382     self.assertTrue( abs(d2-d0-integrate(dg0[1]*p)) < 1e-2 * abs(d2-d0) )
383    
384 gross 4688
385     if __name__ == "__main__":
386     suite = unittest.TestSuite()
387 gross 4979 #suite.addTest(unittest.makeSuite(TestAcousticInversion))
388     suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))
389 gross 4688 s=unittest.TextTestRunner(verbosity=2).run(suite)
390     if not s.wasSuccessful(): sys.exit(1)
391    

  ViewVC Help
Powered by ViewVC 1.1.26