/[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 5556 by gross, Wed Mar 25 00:11:28 2015 UTC revision 5558 by caltinay, Wed Mar 25 02:13:05 2015 UTC
# Line 128  class TestAcousticInversion(unittest.Tes Line 128  class TestAcousticInversion(unittest.Tes
128          self.assertTrue(isinstance(dg, Data))          self.assertTrue(isinstance(dg, Data))
129          self.assertTrue(dg.getShape()==(2,))          self.assertTrue(dg.getShape()==(2,))
130          self.assertTrue(dg.getFunctionSpace()==Solution(domain))          self.assertTrue(dg.getFunctionSpace()==Solution(domain))
131          self.assertTrue(Lsup(dg) < 1e-10)          self.assertLess(Lsup(dg), 1e-10)
132    
133          # this shuld be zero' too          # this shuld be zero' too
134          sigma_comps=[2*sigma.real, sigma.imag/2.]          sigma_comps=[2*sigma.real, sigma.imag/2.]
135          args=acw.getArguments(sigma_comps)          args=acw.getArguments(sigma_comps)
136          d=acw.getDefect(sigma_comps, *args)          d=acw.getDefect(sigma_comps, *args)
137          self.assertTrue(isinstance(d, float))          self.assertTrue(isinstance(d, float))
138          self.assertTrue(abs(d)< 1e-10)          self.assertLess(abs(d), 1e-10)
139    
140          dg=acw.getGradient(sigma_comps, *args)          dg=acw.getGradient(sigma_comps, *args)
141          self.assertTrue(isinstance(dg, Data))          self.assertTrue(isinstance(dg, Data))
142          self.assertTrue(dg.getShape()==(2,))          self.assertTrue(dg.getShape()==(2,))
143          self.assertTrue(dg.getFunctionSpace()==Solution(domain))          self.assertTrue(dg.getFunctionSpace()==Solution(domain))
144          self.assertTrue(Lsup(dg) < 1e-10)          self.assertLess(Lsup(dg), 1e-10)
145    
146          # this shouldn't be zero:          # this shouldn't be zero:
147          sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])          sigma0=[2*sigma.real, 10*a.imag]*(27*Function(domain).getX()[0]-Function(domain).getX()[1])
# Line 166  class TestAcousticInversion(unittest.Tes Line 166  class TestAcousticInversion(unittest.Tes
166          sigma1=sigma0+p*[1,0]          sigma1=sigma0+p*[1,0]
167          args=acw.getArguments(sigma1)          args=acw.getArguments(sigma1)
168          d1=acw.getDefect(sigma1, *args)          d1=acw.getDefect(sigma1, *args)
169          self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
170    
171          sigma2=sigma0+p*[0,1]          sigma2=sigma0+p*[0,1]
172          args=acw.getArguments(sigma2)          args=acw.getArguments(sigma2)
173          d2=acw.getDefect(sigma2, *args)          d2=acw.getDefect(sigma2, *args)
174          self.assertTrue( abs(d2-d0-integrate(dg0[1]*p))  < 1e-2  * abs(d2-d0) )          self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
175    
176      def test_numeric2DnoscaleF(self):      def test_numeric2DnoscaleF(self):
177          domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])          domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
# Line 193  class TestAcousticInversion(unittest.Tes Line 193  class TestAcousticInversion(unittest.Tes
193          args=acw.getArguments(sigma_comps)          args=acw.getArguments(sigma_comps)
194          d=acw.getDefect(sigma_comps, *args)          d=acw.getDefect(sigma_comps, *args)
195          self.assertTrue(isinstance(d, float))          self.assertTrue(isinstance(d, float))
196          self.assertTrue(Lsup(d) < 1e-10)          self.assertLess(Lsup(d), 1e-10)
197          #self.assertTrue(d >= 0)          #self.assertTrue(d >= 0)
         #self.assertTrue(d < 1e-10)  
198    
199          dg=acw.getGradient(sigma_comps, *args)          dg=acw.getGradient(sigma_comps, *args)
200    
201          self.assertTrue(isinstance(dg, Data))          self.assertTrue(isinstance(dg, Data))
202          self.assertTrue(dg.getShape()==(2,))          self.assertTrue(dg.getShape()==(2,))
203          self.assertTrue(dg.getFunctionSpace()==Solution(domain))          self.assertTrue(dg.getFunctionSpace()==Solution(domain))
204          self.assertTrue(Lsup(dg) < 5e-10)          self.assertLess(Lsup(dg), 5e-10)
205          # this shouldn't be zero:          # this shouldn't be zero:
206          sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )          sigma0=Data([2*sigma.real, sigma.imag/2], Function(domain) )
207          args=acw.getArguments(sigma0)          args=acw.getArguments(sigma0)
# Line 225  class TestAcousticInversion(unittest.Tes Line 224  class TestAcousticInversion(unittest.Tes
224          args=acw.getArguments(sigma1)          args=acw.getArguments(sigma1)
225          d1=acw.getDefect(sigma1, *args)          d1=acw.getDefect(sigma1, *args)
226    
227          self.assertTrue( abs( d1-d0-integrate(dg0[0]*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0[0]*p) ), 1e-2*abs(d1-d0) )
228    
229          sigma2=sigma0+p*[0,1]          sigma2=sigma0+p*[0,1]
230          args=acw.getArguments(sigma2)          args=acw.getArguments(sigma2)
231          d2=acw.getDefect(sigma2, *args)          d2=acw.getDefect(sigma2, *args)
232          self.assertTrue( abs(d2-d0-integrate(dg0[1]*p))  < 1e-2  * abs(d2-d0) )          self.assertLess( abs(d2-d0-integrate(dg0[1]*p)), 1e-2*abs(d2-d0) )
233    
234    
235  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
# Line 250  class TestSubsidence(unittest.TestCase): Line 249  class TestSubsidence(unittest.TestCase):
249          P0=10.          P0=10.
250          args0=acw.getArguments(P0)          args0=acw.getArguments(P0)
251          u=args0[0]          u=args0[0]
252          self.assertTrue(Lsup(u[0]) < 1.e-8)          self.assertLess(Lsup(u[0]), 1.e-8)
253          self.assertTrue(Lsup(u[1]) < 1.e-8)          self.assertLess(Lsup(u[1]), 1.e-8)
254          self.assertTrue(Lsup(u[2]-2.5*domain.getX()[2]) < 1.e-8)          self.assertLess(Lsup(u[2]-2.5*domain.getX()[2]), 1.e-8)
255    
256          dd=acw.getDefect(P0, *args0)          dd=acw.getDefect(P0, *args0)
257    
# Line 282  class TestSubsidence(unittest.TestCase): Line 281  class TestSubsidence(unittest.TestCase):
281          args1=acw.getArguments(P1)          args1=acw.getArguments(P1)
282          d1=acw.getDefect(P1, *args1)          d1=acw.getDefect(P1, *args1)
283          ref=abs((d1-d0)/INC)          ref=abs((d1-d0)/INC)
284          self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dP)) < ref * 1.e-5)          self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dP)), ref * 1.e-5)
285    
286          dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)          dP=exp(-(length(x-[0.3,0.3,0.5])/0.06)**2)
287          P2=P0-INC*dP          P2=P0-INC*dP
288          args2=acw.getArguments(P2)          args2=acw.getArguments(P2)
289          d2=acw.getDefect(P2, *args2)          d2=acw.getDefect(P2, *args2)
290          ref=abs((d2-d0)/INC)          ref=abs((d2-d0)/INC)
291          self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)          self.assertLess(abs((d2-d0)/INC+integrate(grad_d* dP)), ref * 1.e-5)
292    
293  @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")  @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
294  class TestDCResistivity(unittest.TestCase):  class TestDCResistivity(unittest.TestCase):
# Line 325  class TestDCResistivity(unittest.TestCas Line 324  class TestDCResistivity(unittest.TestCas
324    
325          acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)          acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)
326    
327          self.assertTrue(Lsup(phiPrimary-acw.getPrimaryPotential()) < 1.e-10 * Lsup(acw.getPrimaryPotential()))          self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))
328    
329          SIGMA=10. # matches current          SIGMA=10. # matches current
330          args0=acw.getArguments(SIGMA)          args0=acw.getArguments(SIGMA)
# Line 334  class TestDCResistivity(unittest.TestCas Line 333  class TestDCResistivity(unittest.TestCas
333    
334          # true secondary potential          # true secondary potential
335          pps=pp-phiPrimary          pps=pp-phiPrimary
336          self.assertTrue(Lsup(p-pps) < 1.e-6 * Lsup(pps))          self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))
   
337    
338          # test return values at electrodes:          # test return values at electrodes:
339          self.assertTrue(abs(u[0]-uu[0]*uuscale) < 1.e-6 * abs(uu[0]*uuscale))          self.assertLess(abs(u[0]-uu[0]*uuscale), 1.e-6 * abs(uu[0]*uuscale))
340          self.assertTrue(abs(u[1]-uu[1]*uuscale) < 1.e-6 * abs(uu[1]*uuscale))          self.assertLess(abs(u[1]-uu[1]*uuscale), 1.e-6 * abs(uu[1]*uuscale))
341    
342          # this sould be zero          # this sould be zero
343          dd=acw.getDefect(SIGMA, *args0)          dd=acw.getDefect(SIGMA, *args0)
# Line 378  class TestDCResistivity(unittest.TestCas Line 376  class TestDCResistivity(unittest.TestCas
376          args1=acw.getArguments(SIGMA1)          args1=acw.getArguments(SIGMA1)
377          d1=acw.getDefect(SIGMA1, *args1)          d1=acw.getDefect(SIGMA1, *args1)
378          ref=abs((d1-d0)/INC)          ref=abs((d1-d0)/INC)
379          self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)          self.assertLess(abs((d1-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
380    
381          dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)          dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
382          SIGMA2=SIGMA0+INC*dS          SIGMA2=SIGMA0+INC*dS
383          args2=acw.getArguments(SIGMA2)          args2=acw.getArguments(SIGMA2)
384          d2=acw.getDefect(SIGMA2, *args2)          d2=acw.getDefect(SIGMA2, *args2)
385          ref=abs((d2-d0)/INC)          ref=abs((d2-d0)/INC)
386          self.assertTrue(abs((d2-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)          self.assertLess(abs((d2-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
387    
388          dS=-1          dS=-1
389          SIGMA3=SIGMA0+INC*dS          SIGMA3=SIGMA0+INC*dS
390          args3=acw.getArguments(SIGMA3)          args3=acw.getArguments(SIGMA3)
391          d3=acw.getDefect(SIGMA3, *args3)          d3=acw.getDefect(SIGMA3, *args3)
392          ref=abs((d3-d0)/INC)          ref=abs((d3-d0)/INC)
393          self.assertTrue(abs((d3-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)          self.assertLess(abs((d3-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
394    
395          dS=1          dS=1
396          SIGMA4=SIGMA0+INC*dS          SIGMA4=SIGMA0+INC*dS
397          args4=acw.getArguments(SIGMA4)          args4=acw.getArguments(SIGMA4)
398          d4=acw.getDefect(SIGMA4, *args4)          d4=acw.getDefect(SIGMA4, *args4)
399          ref=abs((d4-d0)/INC)          ref=abs((d4-d0)/INC)
400          self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)          self.assertLess(abs((d4-d0)/INC-integrate(grad_d* dS)), ref * 1.e-3)
401    
402  class TestIsostaticPressure(unittest.TestCase):  class TestIsostaticPressure(unittest.TestCase):
403      @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")      @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
# Line 411  class TestIsostaticPressure(unittest.Tes Line 409  class TestIsostaticPressure(unittest.Tes
409          rho=Scalar(100, Function(domain))          rho=Scalar(100, Function(domain))
410          p0=ps.getPressure(g, rho)          p0=ps.getPressure(g, rho)
411          p_ref=-(1.-domain.getX()[2])*981.          p_ref=-(1.-domain.getX()[2])*981.
412          self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))          self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
413    
414          g=Vector([0,0,-10], Function(domain))          g=Vector([0,0,-10], Function(domain))
415          rho=Scalar(0, Function(domain))          rho=Scalar(0, Function(domain))
416          p0=ps.getPressure(g, rho)          p0=ps.getPressure(g, rho)
417          p_ref=-(1.-domain.getX()[2])*26700          p_ref=-(1.-domain.getX()[2])*26700
418          self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))          self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
419    
420          g=Vector([0,0,-10], Function(domain))          g=Vector([0,0,-10], Function(domain))
421          rho=Scalar(100, Function(domain))          rho=Scalar(100, Function(domain))
422          p0=ps.getPressure(g, rho)          p0=ps.getPressure(g, rho)
423          p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)          p_ref=-(1.-domain.getX()[2])*(981.+26700+1000)
424          self.assertTrue(Lsup(p0-p_ref) < 1e-6 * Lsup(p_ref))          self.assertLess(Lsup(p0-p_ref), 1e-6 * Lsup(p_ref))
425    
426  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
427  class TestMT2DModelTEMode(unittest.TestCase):  class TestMT2DModelTEMode(unittest.TestCase):
# Line 436  class TestMT2DModelTEMode(unittest.TestC Line 434  class TestMT2DModelTEMode(unittest.TestC
434          w0=1.          w0=1.
435          Ex0=1.          Ex0=1.
436          # now we do a real one          # now we do a real one
437          acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, Ex_top=Ex0)          model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, w0=w0, Ex_top=Ex0)
438          self.assertEqual(acw.getDomain(),  domain)          self.assertEqual(model.getDomain(),  domain)
439          pde=acw.setUpPDE()          pde=model.setUpPDE()
440          self.assertIsInstance(pde, LinearPDE)          self.assertIsInstance(pde, LinearPDE)
441          self.assertEqual(pde.getNumEquations(), 2)          self.assertEqual(pde.getNumEquations(), 2)
442          self.assertEqual(pde.getNumSolutions(), 2)          self.assertEqual(pde.getNumSolutions(), 2)
443          self.assertEqual(pde.getDomain(),  domain)          self.assertEqual(pde.getDomain(),  domain)
444    
445          # other things that should work          # other things that should work
446          acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Ex_top=complex(4.5,6) )          model=MT2DModelTEMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Ex_top=complex(4.5,6) )
447    
448          # these shouldn't work          # these shouldn't work
449          self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )          self.assertRaises(ValueError, MT2DModelTEMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Ex_top=complex(4.5,6) )
# Line 477  class TestMT2DModelTEMode(unittest.TestC Line 475  class TestMT2DModelTEMode(unittest.TestC
475          Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))          Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
476          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real
477    
478          acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9,  directSolver=True)          model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9,  directSolver=True)
479    
480          args=acw.getArguments(SIGMA)          args=model.getArguments(SIGMA)
481          Ex=args[0]          Ex=args[0]
482          Exz=args[1]          Exz=args[1]
483          self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))          self.assertTrue(Lsup(Ex[0]-Ex0_ex) <= 1e-4 * Lsup(Ex0_ex))
# Line 487  class TestMT2DModelTEMode(unittest.TestC Line 485  class TestMT2DModelTEMode(unittest.TestC
485          self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))          self.assertTrue(Lsup(Exz[0]-Ex0_ex_z) <= 1e-2 * Lsup(Ex0_ex_z))
486          self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))          self.assertTrue(Lsup(Exz[1]-Ex1_ex_z) <= 1e-2 * Lsup(Ex1_ex_z))
487    
488          argsr=acw.getArguments(0.)          argsr=model.getArguments(0.)
489          ref=acw.getDefect(0., *argsr)          ref=model.getDefect(0., *argsr)
490    
491          # this should be almost zero:          # this should be almost zero:
492          args=acw.getArguments(SIGMA)          args=model.getArguments(SIGMA)
493          d=acw.getDefect(SIGMA, *args)          d=model.getDefect(SIGMA, *args)
494          self.assertTrue( d > 0.)          self.assertTrue( d > 0.)
495          self.assertTrue( ref > 0.)          self.assertTrue( ref > 0.)
496          self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)          self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
# Line 503  class TestMT2DModelTEMode(unittest.TestC Line 501  class TestMT2DModelTEMode(unittest.TestC
501          Ex0_ex_z=-sin(k.imag*z)*k.imag*(exp(k.real*z)-exp(-k.real*z))+cos(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))*k.real          Ex0_ex_z=-sin(k.imag*z)*k.imag*(exp(k.real*z)-exp(-k.real*z))+cos(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))*k.real
502          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real
503          # and this should be zero          # and this should be zero
504          d0=acw.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])          d0=model.getDefect(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
505          self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)          self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)
506    
507          # and this too          # and this too
508          dg=acw.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])          dg=model.getGradient(SIGMA, Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], Ex0_ex_z*[1.,0]+ Ex1_ex_z*[0,1.])
509          self.assertTrue(isinstance(dg, Data))          self.assertTrue(isinstance(dg, Data))
510          self.assertTrue(dg.getShape()==())          self.assertTrue(dg.getShape()==())
511          self.assertTrue(Lsup(dg) < 1e-10)          self.assertLess(Lsup(dg), 1e-10)
512    
513      def test_Differential(self):      def test_Differential(self):
514          INC=0.001          INC=0.001
# Line 537  class TestMT2DModelTEMode(unittest.TestC Line 535  class TestMT2DModelTEMode(unittest.TestC
535          Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))          Ex1_ex=sin(k.imag*z)*(exp(k.real*z)+exp(-k.real*z))
536          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real          Ex1_ex_z=cos(k.imag*z)*k.imag*(exp(k.real*z)+exp(-k.real*z))+sin(k.imag*z)*(exp(k.real*z)-exp(-k.real*z))*k.real
537    
538          acw=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9,  directSolver=True)          model=MT2DModelTEMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtTop=True, Ex_top=Ex0_ex*[1.,0]+ Ex1_ex*[0,1.], tol=1e-9,  directSolver=True)
539    
540          # this is the base line:          # this is the base line:
541          xx=domain.getX()[0]          xx=domain.getX()[0]
542          SIGMA0=3.*(xx+0.3)          SIGMA0=3.*(xx+0.3)
543          args0=acw.getArguments(SIGMA0)          args0=model.getArguments(SIGMA0)
544          d0=acw.getDefect(SIGMA0, *args0)          d0=model.getDefect(SIGMA0, *args0)
545          dg0=acw.getGradient(SIGMA0, *args0)          dg0=model.getGradient(SIGMA0, *args0)
546          self.assertTrue(isinstance(dg0, Data))          self.assertTrue(isinstance(dg0, Data))
547          self.assertTrue(dg0.getShape()==())          self.assertTrue(dg0.getShape()==())
548    
# Line 553  class TestMT2DModelTEMode(unittest.TestC Line 551  class TestMT2DModelTEMode(unittest.TestC
551          # test 1          # test 1
552          p=INC          p=INC
553          SIGMA1=SIGMA0+p          SIGMA1=SIGMA0+p
554          args1=acw.getArguments(SIGMA1)          args1=model.getArguments(SIGMA1)
555          d1=acw.getDefect(SIGMA1, *args1)          d1=model.getDefect(SIGMA1, *args1)
556          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
557    
558          # test 2          # test 2
559          p=exp(-length(X-(0.2,0.2))**2/10)*INC          p=exp(-length(X-(0.2,0.2))**2/10)*INC
560          SIGMA1=SIGMA0+p          SIGMA1=SIGMA0+p
561          args1=acw.getArguments(SIGMA1)          args1=model.getArguments(SIGMA1)
562          d1=acw.getDefect(SIGMA1, *args1)          d1=model.getDefect(SIGMA1, *args1)
563          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
564    
565          # test 3          # test 3
566          p=sin(length(X)*3*3.14)*INC          p=sin(length(X)*3*3.14)*INC
567          SIGMA1=SIGMA0+p          SIGMA1=SIGMA0+p
568          args1=acw.getArguments(SIGMA1)          args1=model.getArguments(SIGMA1)
569          d1=acw.getDefect(SIGMA1, *args1)          d1=model.getDefect(SIGMA1, *args1)
570          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
571    
572    
573  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")  @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
# Line 583  class TestMT2DModelTMMode(unittest.TestC Line 581  class TestMT2DModelTMMode(unittest.TestC
581          w0=1.          w0=1.
582          Hx0=1.          Hx0=1.
583          # now we do a real one          # now we do a real one
584          acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)          model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, w0=w0, Hx_bottom=Hx0)
585          self.assertEqual(acw.getDomain(),  domain)          self.assertEqual(model.getDomain(),  domain)
586          pde=acw.setUpPDE()          pde=model.setUpPDE()
587          self.assertIsInstance(pde, LinearPDE)          self.assertIsInstance(pde, LinearPDE)
588          self.assertEqual(pde.getNumEquations(), 2)          self.assertEqual(pde.getNumEquations(), 2)
589          self.assertEqual(pde.getNumSolutions(), 2)          self.assertEqual(pde.getNumSolutions(), 2)
590          self.assertEqual(pde.getDomain(),  domain)          self.assertEqual(pde.getDomain(),  domain)
591    
592          # other things that should work          # other things that should work
593          acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )          model=MT2DModelTMMode(domain, omega, x, Z_XY, eta=None, w0=[2.,3.], Hx_bottom=complex(4.5,6) )
594    
595          # these shouldn't work          # these shouldn't work
596          self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )          self.assertRaises(ValueError, MT2DModelTMMode, domain, omega, x, [3.], eta=[1.,1.], w0=[2.,3.], Hx_bottom=complex(4.5,6) )
# Line 625  class TestMT2DModelTMMode(unittest.TestC Line 623  class TestMT2DModelTMMode(unittest.TestC
623          Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2          Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
624          Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2          Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2
625    
626          acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, Hx_bottom=Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], tol=1e-9,  directSolver=True)          model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, fixAtBottom=True, Hx_bottom=Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], tol=1e-9,  directSolver=True)
627    
628          args=acw.getArguments(RHO)          args=model.getArguments(RHO)
629          Hx=args[0]          Hx=args[0]
630          Hxz=args[1]          g_Hx=args[1]
631          self.assertTrue(Lsup(Hx[0]-Hx0_ex) <= 1e-4 * Lsup(Hx0_ex))          Hxz=g_Hx[:,1]
632          self.assertTrue(Lsup(Hx[1]-Hx1_ex) <= 1e-4 * Lsup(Hx1_ex))          self.assertLess(Lsup(Hx[0]-Hx0_ex), 1e-4 * Lsup(Hx0_ex))
633          self.assertTrue(Lsup(Hxz[0]-Hx0_ex_z) <= 1e-2 * Lsup(Hx0_ex_z))          self.assertLess(Lsup(Hx[1]-Hx1_ex), 1e-4 * Lsup(Hx1_ex))
634          self.assertTrue(Lsup(Hxz[1]-Hx1_ex_z) <= 1e-2 * Lsup(Hx1_ex_z))          self.assertLess(Lsup(Hxz[0]-Hx0_ex_z), 1e-2 * Lsup(Hx0_ex_z))
635            self.assertLess(Lsup(Hxz[1]-Hx1_ex_z), 1e-2 * Lsup(Hx1_ex_z))
636    
637          argsr=acw.getArguments(1.)          argsr=model.getArguments(1.)
638          ref=acw.getDefect(1., *argsr)          ref=model.getDefect(1., *argsr)
639    
640          # this should be almost zero:          # this should be almost zero:
641          args=acw.getArguments(RHO)          args=model.getArguments(RHO)
642          d=acw.getDefect(RHO, *args)          d=model.getDefect(RHO, *args)
643          self.assertTrue( d > 0.)          self.assertTrue( d > 0.)
644          self.assertTrue( ref > 0.)          self.assertTrue( ref > 0.)
645          self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)          self.assertTrue( d <= 3e-3 * ref ) # d should be zero (some sort of)
# Line 650  class TestMT2DModelTMMode(unittest.TestC Line 649  class TestMT2DModelTMMode(unittest.TestC
649          Hx0_ex_z=(-sin(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))+exp(-k.real*(z-L)))+cos(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))*k.real)/2          Hx0_ex_z=(-sin(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))+exp(-k.real*(z-L)))+cos(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))*k.real)/2
650          Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2          Hx1_ex=sin(k.imag*(z-L))*(exp(k.real*(z-L))-exp(-k.real*(z-L)))/2
651          Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2          Hx1_ex_z=(cos(k.imag*(z-L))*k.imag*(exp(k.real*(z-L))-exp(-k.real*(z-L)))+sin(k.imag*(z-L))*(exp(k.real*(z-L))+exp(-k.real*(z-L)))*k.real)/2
652            g_Hx = Data(0, (2,2), Hx0_ex_z.getFunctionSpace())
653            g_Hx[0,1] = Hx0_ex_z
654            g_Hx[1,1] = Hx1_ex_z
655          # and this should be zero          # and this should be zero
656          d0=acw.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])          d0=model.getDefect(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], g_Hx)
657          self.assertTrue( d0 <= 1e-8 * ref ) # d should be zero (some sort of)          self.assertLess( d0, 1e-8 * ref ) # d should be zero (some sort of)
658    
659          # and this too          # and this too
660          dg=acw.getGradient(RHO, Hx0_ex*[1.,0]+ Hx1_ex*[0,1.], Hx0_ex_z*[1.,0]+ Hx1_ex_z*[0,1.])          dg=model.getGradient(RHO, Hx0_ex*[1.,0]+Hx1_ex*[0,1.], g_Hx)
661          self.assertTrue(isinstance(dg, Data))          self.assertTrue(isinstance(dg, Data))
662          self.assertTrue(dg.getShape()==())          self.assertTrue(dg.getShape()==())
663          self.assertTrue(Lsup(dg) < 1e-10)          self.assertLess(Lsup(dg), 1e-10)
664    
665      def test_Differential(self):      def test_Differential(self):
666          INC=0.001          INC=0.001
# Line 680  class TestMT2DModelTMMode(unittest.TestC Line 682  class TestMT2DModelTMMode(unittest.TestC
682          x=[ [X1,Z0], [X2,Z0] ]          x=[ [X1,Z0], [X2,Z0] ]
683          eta=None          eta=None
684    
685          acw=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9,  directSolver=True)          model=MT2DModelTMMode(domain, omega, x, Z_XY, eta, mu=mu0, tol=1e-9,  directSolver=True)
686    
687          # this is the base line:          # this is the base line:
688          xx=domain.getX()[0]          xx=domain.getX()[0]
689          RHO0=3.*(xx+0.3)          RHO0=3.*(xx+0.3)
690          args0=acw.getArguments(RHO0)          args0=model.getArguments(RHO0)
691          d0=acw.getDefect(RHO0, *args0)          d0=model.getDefect(RHO0, *args0)
692          dg0=acw.getGradient(RHO0, *args0)          dg0=model.getGradient(RHO0, *args0)
693          self.assertTrue(isinstance(dg0, Data))          self.assertTrue(isinstance(dg0, Data))
694          self.assertTrue(dg0.getShape()==())          self.assertTrue(dg0.getShape()==())
695    
# Line 696  class TestMT2DModelTMMode(unittest.TestC Line 698  class TestMT2DModelTMMode(unittest.TestC
698          # test 1          # test 1
699          p=INC          p=INC
700          RHO1=RHO0+p          RHO1=RHO0+p
701          args1=acw.getArguments(RHO1)          args1=model.getArguments(RHO1)
702          d1=acw.getDefect(RHO1, *args1)          d1=model.getDefect(RHO1, *args1)
703          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
704    
705          # test 2          # test 2
706          p=exp(-length(X-(0.2,0.2))**2/10)*INC          p=exp(-length(X-(0.2,0.2))**2/10)*INC
707          RHO1=RHO0+p          RHO1=RHO0+p
708          args1=acw.getArguments(RHO1)          args1=model.getArguments(RHO1)
709          d1=acw.getDefect(RHO1, *args1)          d1=model.getDefect(RHO1, *args1)
710          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
711    
712          # test 3          # test 3
713          p=sin(length(X)*3*3.14)*INC          p=sin(length(X)*3*3.14)*INC
714          RHO1=RHO0+p          RHO1=RHO0+p
715          args1=acw.getArguments(RHO1)          args1=model.getArguments(RHO1)
716          d1=acw.getDefect(RHO1, *args1)          d1=model.getDefect(RHO1, *args1)
717          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertLess( abs( d1-d0-integrate(dg0*p) ), 1e-2*abs(d1-d0) )
718    
719    
720  if __name__ == '__main__':  if __name__ == '__main__':
721      run_tests(__name__, exit_on_failure=True)      run_tests(__name__, exit_on_failure=True)
722    

Legend:
Removed from v.5556  
changed lines
  Added in v.5558

  ViewVC Help
Powered by ViewVC 1.1.26