/[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 5262 by gross, Thu Nov 13 08:11:40 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1 
1    from __future__ import print_function
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2014 by University of Queensland  # Copyright (c) 2003-2014 by University of Queensland
# Line 36  from esys.escript.linearPDEs import Line Line 36  from esys.escript.linearPDEs import Line
36  from esys.escript import getEscriptParamInt  from esys.escript import getEscriptParamInt
37  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
38    
39    try:
40        from esys.ripley import Rectangle, Brick as ripRectangle, ripBrick
41        HAVE_RIPLEY = True
42    except ImportError as e:
43        HAVE_RIPLEY = False
44    
45    try:
46        from esys.finley import Rectangle as ripRectangle
47        HAVE_FINLEY = True
48    except ImportError as e:
49        HAVE_FINLEY = False
50    
51  mpisize = getMPISizeWorld()  mpisize = getMPISizeWorld()
52  # this is mainly to avoid warning messages  # this is mainly to avoid warning messages
53  logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)  logging.basicConfig(format='%(name)s: %(message)s', level=logging.INFO)
# Line 53  except KeyError: Line 65  except KeyError:
65            
66  have_direct=getEscriptParamInt("PASO_DIRECT")    have_direct=getEscriptParamInt("PASO_DIRECT")  
67    
68    @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
69  @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")  @unittest.skipIf(mpisize>1 or have_direct!=1, "more than 1 MPI rank or missing direct solver")
70  class TestAcousticInversion(unittest.TestCase):  class TestAcousticInversion(unittest.TestCase):
71      def test_API(self):      def test_API(self):
72          from esys.ripley import Rectangle  
73          domain=Rectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])          domain=ripRectangle(20,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
74          omega=2.          omega=2.
75                    
76                    
# Line 84  class TestAcousticInversion(unittest.Tes Line 96  class TestAcousticInversion(unittest.Tes
96                    
97      def test_numeric2DscaleF(self):      def test_numeric2DscaleF(self):
98                    
99          from esys.ripley import Rectangle          domain=ripRectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])
         domain=Rectangle(100,100, diracPoints=[(0.5,1.)], diracTags=['sss'])  
100          omega=2.          omega=2.
101                    
102          # test solution is u = a * z where a is complex          # test solution is u = a * z where a is complex
# Line 169  class TestAcousticInversion(unittest.Tes Line 180  class TestAcousticInversion(unittest.Tes
180            
181      def test_numeric2DnoscaleF(self):      def test_numeric2DnoscaleF(self):
182                    
183          from esys.ripley import Rectangle          domain=ripRectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])
         domain=Rectangle(10,20, diracPoints=[(0.5,1.)], diracTags=['sss'])  
184          omega=1.5          omega=1.5
185                    
186          # test solution is u = a * z where a is complex          # test solution is u = a * z where a is complex
# Line 229  class TestAcousticInversion(unittest.Tes Line 239  class TestAcousticInversion(unittest.Tes
239          d2=acw.getDefect(sigma2, *args)          d2=acw.getDefect(sigma2, *args)
240          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) )
241    
242    @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
243  class TestMT2DModelTEMode(unittest.TestCase):  class TestMT2DModelTEMode(unittest.TestCase):
244      def test_API(self):      def test_API(self):
245          from esys.ripley import Rectangle          domain=ripRectangle(19, 19, d1=mpisize)
         domain=Rectangle(19, 19, d1=mpisize)  
246          omega=2.          omega=2.
247          x=[ [0.2,0.5], [0.3,0.5] ]          x=[ [0.2,0.5], [0.3,0.5] ]
248          Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]          Z_XY=[ complex(1.2,1.5), complex(1.3,2.5) ]
# Line 262  class TestMT2DModelTEMode(unittest.TestC Line 272  class TestMT2DModelTEMode(unittest.TestC
272          SIGMA=15.          SIGMA=15.
273          k=cmath.sqrt(1j*omega*mu0*SIGMA)  # Ex=exp(k*z)          k=cmath.sqrt(1j*omega*mu0*SIGMA)  # Ex=exp(k*z)
274    
275          from esys.ripley import Rectangle          domain=ripRectangle(199,199, d1=mpisize)
         domain=Rectangle(199,199, d1=mpisize)  
276    
277                    
278          IMP=cmath.sqrt(1j*omega*mu0/SIGMA)          IMP=cmath.sqrt(1j*omega*mu0/SIGMA)
# Line 317  class TestMT2DModelTEMode(unittest.TestC Line 326  class TestMT2DModelTEMode(unittest.TestC
326          SIGMA=15.          SIGMA=15.
327          k=cmath.sqrt(1j*omega*mu0*SIGMA)  # Ex=exp(k*z)          k=cmath.sqrt(1j*omega*mu0*SIGMA)  # Ex=exp(k*z)
328    
329          from esys.ripley import Rectangle          domain=ripRectangle(99,99, d1=mpisize)
         domain=Rectangle(99,99, d1=mpisize)  
330    
331                    
332          IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)          IMP=-cmath.sqrt(1j*omega*mu0/SIGMA)
# Line 366  class TestMT2DModelTEMode(unittest.TestC Line 374  class TestMT2DModelTEMode(unittest.TestC
374          d1=acw.getDefect(SIGMA1, *args1)          d1=acw.getDefect(SIGMA1, *args1)
375          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )          self.assertTrue( abs( d1-d0-integrate(dg0*p) ) < 1e-2  * abs(d1-d0) )
376    
377    @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
378  class TestSubsidence(unittest.TestCase):  class TestSubsidence(unittest.TestCase):
379      def test_PDE(self):      def test_PDE(self):
380                    
381          lam=2.          lam=2.
382          mu=1.          mu=1.
383    
384          from esys.ripley import Brick          domain=ripBrick(20,20,19, d2=mpisize)
         domain=Brick(20,20,19, d2=mpisize)  
385                    
386          xb=FunctionOnBoundary(domain).getX()          xb=FunctionOnBoundary(domain).getX()
387          m=whereZero(xb[2]-1)          m=whereZero(xb[2]-1)
# Line 398  class TestSubsidence(unittest.TestCase): Line 406  class TestSubsidence(unittest.TestCase):
406          mu=1.          mu=1.
407                    
408          INC=0.01          INC=0.01
409          from esys.ripley import Brick          domain=ripBrick(20,20,20*mpisize-1 , d2=mpisize)
         domain=Brick(20,20,20*mpisize-1 , d2=mpisize)  
410                    
411          xb=FunctionOnBoundary(domain).getX()          xb=FunctionOnBoundary(domain).getX()
412          m=whereZero(xb[2]-1)          m=whereZero(xb[2]-1)
# Line 429  class TestSubsidence(unittest.TestCase): Line 436  class TestSubsidence(unittest.TestCase):
436          ref=abs((d2-d0)/INC)          ref=abs((d2-d0)/INC)
437          self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)          self.assertTrue(abs((d2-d0)/INC+integrate(grad_d* dP)) < ref * 1.e-5)
438    
439    @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
440  class TestDCResistivity(unittest.TestCase):  class TestDCResistivity(unittest.TestCase):
441    
442      def test_PDE2D(self):      def test_PDE2D(self):
# Line 437  class TestDCResistivity(unittest.TestCas Line 445  class TestDCResistivity(unittest.TestCas
445    
446          sigma0=1.          sigma0=1.
447          electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]          electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
448          from esys.finley import Rectangle          domain=finRectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
         domain=Rectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )  
449          loc=Locator(domain,electrodes[2:])          loc=Locator(domain,electrodes[2:])
450    
451          # this creates some reference Data:          # this creates some reference Data:
# Line 493  class TestDCResistivity(unittest.TestCas Line 500  class TestDCResistivity(unittest.TestCas
500          sigma0=1.          sigma0=1.
501          dx_tests=0.1          dx_tests=0.1
502          electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]          electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
503          from esys.finley import Rectangle          domain=finRectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
         domain=Rectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )  
504          loc=Locator(domain,electrodes[2:])          loc=Locator(domain,electrodes[2:])
505    
506          # arguments for DcRes          # arguments for DcRes
# Line 545  class TestDCResistivity(unittest.TestCas Line 551  class TestDCResistivity(unittest.TestCas
551          self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)          self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
552    
553  class TestIsostaticPressure(unittest.TestCase):  class TestIsostaticPressure(unittest.TestCase):
554        @unittest.skipIf(not HAVE_RIPLEY, "Ripley module not available")
555      def test_all(self):      def test_all(self):
556          from esys.ripley import Brick          domain=ripBrick(50,50,20*mpisize-1, d2=mpisize)
         domain=Brick(50,50,20*mpisize-1, d2=mpisize)  
557    
558          ps=IsostaticPressure(domain, level0=1., coordinates=None)          ps=IsostaticPressure(domain, level0=1., coordinates=None)
559            

Legend:
Removed from v.5262  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26