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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4938 by jfenwick, Wed May 14 01:13:23 2014 UTC revision 4979 by gross, Sun Jun 1 22:36:16 2014 UTC
# Line 26  import esys.escriptcore.utestselect as u Line 26  import esys.escriptcore.utestselect as u
26  import numpy as np  import numpy as np
27  import os  import os
28  import sys  import sys
29    import cmath
30  from esys.downunder import *  from esys.downunder import *
31  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
32  from esys.escript import *  from esys.escript import *
# Line 226  class TestAcousticInversion(unittest.Tes Line 227  class TestAcousticInversion(unittest.Tes
227          d2=acw.getDefect(sigma2, *args)          d2=acw.getDefect(sigma2, *args)
228          self.assertTrue( abs(d2-d0-integrate(dg0[1]*p))  < 1e-2  * abs(d2-d0) )          self.assertTrue( abs(d2-d0-integrate(dg0[1]*p))  < 1e-2  * abs(d2-d0) )
229    
230    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                                                                        
385  if __name__ == "__main__":  if __name__ == "__main__":
386      suite = unittest.TestSuite()      suite = unittest.TestSuite()
387      suite.addTest(unittest.makeSuite(TestAcousticInversion))      #suite.addTest(unittest.makeSuite(TestAcousticInversion))
388        suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))
389      s=unittest.TextTestRunner(verbosity=2).run(suite)      s=unittest.TextTestRunner(verbosity=2).run(suite)
390      if not s.wasSuccessful(): sys.exit(1)      if not s.wasSuccessful(): sys.exit(1)
391            

Legend:
Removed from v.4938  
changed lines
  Added in v.4979

  ViewVC Help
Powered by ViewVC 1.1.26