/[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 4980 by gross, Sun Jun 1 23:53:11 2014 UTC revision 4981 by gross, Mon Jun 2 00:26:35 2014 UTC
# Line 307  class TestMT2DModelTEMode(unittest.TestC Line 307  class TestMT2DModelTEMode(unittest.TestC
307          self.assertTrue(dg.getShape()==())                        self.assertTrue(dg.getShape()==())              
308          self.assertTrue(Lsup(dg) < 1e-10)          self.assertTrue(Lsup(dg) < 1e-10)
309                    
310      def ttest_Differential(self):      def test_Differential(self):
311            
312          INC=0.001          INC=0.01
313                    
314          omega=2.          omega=2.
315          mu0=0.123          mu0=0.123
# Line 331  class TestMT2DModelTEMode(unittest.TestC Line 331  class TestMT2DModelTEMode(unittest.TestC
331          Ex1_ex_z=(-k.real*sin(-k.imag*z)-k.imag*cos(-k.imag*z)) * exp(-k.real*z)          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.] )          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          # this is the base line:
336            SIGMA0=10.        
337          args0=acw.getArguments(SIGMA0)          args0=acw.getArguments(SIGMA0)
338          d0=acw.getDefect(SIGMA0, *args0)          d0=acw.getDefect(SIGMA0, *args0)
         print "d0 =",d0  
           
         dg=acw.getGradient(SIGMA0, *args0)  
         self.assertTrue(isinstance(dg, Data))  
         self.assertTrue(dg.getShape()==())  
         self.assertTrue(dg.getFunctionSpace()==Solution(domain))  
         print dg                  
         self.assertTrue(Lsup(dg) < 1e-10)  
           
         X=Function(domain).getX()  
         inc=exp(-length(X-(0.2,0.2))**2/0.03)*INC  
           
         SIGMA1=SIGMA0+inc  
         args1=acw.getArguments(SIGMA1)  
         d1=acw.getDefect(SIGMA1, *args1)  
           
           
         print "d0, d1 = ",d0, d1  
   
   
   
   
         1/0  
         # test solution is u = a * z where a is complex  
         a=complex(3.45, 0.56)  
         sigma=complex(1e-3, 0.056)  
           
           
         data=Data([a.real, a.imag], FunctionOnBoundary(domain))  
         mydata=data.copy()  
           
         z=FunctionOnBoundary(domain).getX()[1]  
         w=whereZero(z-1.)  
         # source:  
         F=Data( [1,0],Function(domain))  
         #  
         acw=AcousticWaveForm(domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True)  
         # check rescaled data  
         surv=acw.getSurvey()  
         self.assertAlmostEqual( integrate(length(surv[0])**2 * surv[1]), 1.)  
   
         mydata_scale=sqrt( integrate(w*length(mydata)**2) )  
         self.assertAlmostEqual( acw.getSourceScaling(z*[1, 0.]) , a/mydata_scale )  
         self.assertAlmostEqual( acw.getSourceScaling(mydata) , 1./mydata_scale )  
           
         # this should be zero:  
         sigma_comps=[sigma.real, sigma.imag]  
         args=acw.getArguments(sigma_comps)  
         d=acw.getDefect(sigma_comps, *args)  
         self.assertTrue(isinstance(d, float))  
         self.assertTrue(abs(d) < 1e-10)  
   
         dg=acw.getGradient(sigma_comps, *args)  
         self.assertTrue(isinstance(dg, Data))  
         self.assertTrue(dg.getShape()==(2,))  
         self.assertTrue(dg.getFunctionSpace()==Solution(domain))                  
         self.assertTrue(Lsup(dg) < 1e-10)  
   
         # this shuld be zero' too  
         sigma_comps=[2*sigma.real, sigma.imag/2.]  
         args=acw.getArguments(sigma_comps)  
         d=acw.getDefect(sigma_comps, *args)  
         self.assertTrue(isinstance(d, float))  
         self.assertTrue(abs(d)< 1e-10)  
           
         dg=acw.getGradient(sigma_comps, *args)  
         self.assertTrue(isinstance(dg, Data))  
         self.assertTrue(dg.getShape()==(2,))  
         self.assertTrue(dg.getFunctionSpace()==Solution(domain))                  
         self.assertTrue(Lsup(dg) < 1e-10)  
           
         # this shouldn't be zero:  
         sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])  
         args=acw.getArguments(sigma0)  
         d0=acw.getDefect(sigma0, *args)  
         self.assertTrue(isinstance(d0, float))  
         self.assertTrue(d0 >= 0)  
         self.assertTrue(d0 > 1e-10)  
339                    
340          dg0=acw.getGradient(sigma0, *args)          dg0=acw.getGradient(SIGMA0, *args0)
341          self.assertTrue(isinstance(dg0, Data))          self.assertTrue(isinstance(dg0, Data))
342          self.assertTrue(dg0.getShape()==(2,))          self.assertTrue(dg0.getShape()==())
343          self.assertTrue(dg0.getFunctionSpace()==Solution(domain))          
         self.assertTrue(Lsup(dg0) > 1e-10)  
344                    
         # test the gradient numerrically:  
         h=0.002  
345          X=Function(domain).getX()          X=Function(domain).getX()
         # .. increment:  
         p=h*exp(-(length(X-[0.6,0.6])/10)**2)*Lsup(length(sigma0))  
346    
347                    # test 1
348          sigma1=sigma0+p*[1,0]          p=INC
349          args=acw.getArguments(sigma1)          SIGMA1=SIGMA0+p
350          d1=acw.getDefect(sigma1, *args)          args1=acw.getArguments(SIGMA1)
351          self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2  * abs(d1-d0) )          d1=acw.getDefect(SIGMA1, *args1)
352            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    
361          sigma2=sigma0+p*[0,1]          # test 3
362          args=acw.getArguments(sigma2)          p=sin(length(X)*3*3.14)*INC
363          d2=acw.getDefect(sigma2, *args)          SIGMA1=SIGMA0+p
364          self.assertTrue( abs(d2-d0-integrate(dg0[1]*p))  < 1e-2  * abs(d2-d0) )          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                    
368                                                                        
369  if __name__ == "__main__":  if __name__ == "__main__":
370      suite = unittest.TestSuite()      suite = unittest.TestSuite()
371      #suite.addTest(unittest.makeSuite(TestAcousticInversion))      suite.addTest(unittest.makeSuite(TestAcousticInversion))
372      suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))      suite.addTest(unittest.makeSuite(TestMT2DModelTEMode))
373      s=unittest.TextTestRunner(verbosity=2).run(suite)      s=unittest.TextTestRunner(verbosity=2).run(suite)
374      if not s.wasSuccessful(): sys.exit(1)      if not s.wasSuccessful(): sys.exit(1)

Legend:
Removed from v.4980  
changed lines
  Added in v.4981

  ViewVC Help
Powered by ViewVC 1.1.26