/[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 5235 by caltinay, Wed Oct 29 01:01:44 2014 UTC revision 5262 by gross, Thu Nov 13 08:11:40 2014 UTC
# Line 34  from esys.escript import * Line 34  from esys.escript import *
34  from esys.weipa import saveSilo  from esys.weipa import saveSilo
35  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
36  from esys.escript import getEscriptParamInt  from esys.escript import getEscriptParamInt
37    from esys.escript.pdetools import Locator
38    
39  mpisize = getMPISizeWorld()  mpisize = getMPISizeWorld()
40  # this is mainly to avoid warning messages  # this is mainly to avoid warning messages
# Line 428  class TestSubsidence(unittest.TestCase): Line 429  class TestSubsidence(unittest.TestCase):
429          ref=abs((d2-d0)/INC)          ref=abs((d2-d0)/INC)
430          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)
431    
432    class TestDCResistivity(unittest.TestCase):
433    
434        def test_PDE2D(self):
435            
436            dx_tests=0.1
437    
438            sigma0=1.
439            electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
440            from esys.finley import Rectangle
441            domain=Rectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
442            loc=Locator(domain,electrodes[2:])
443    
444            # this creates some reference Data:
445            x=domain.getX()
446            q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
447            ppde=LinearPDE(domain, numEquations=1)
448            s=Scalar(0.,DiracDeltaFunctions(domain))
449            s.setTaggedValue("sl0" ,1.)
450            s.setTaggedValue("sl1",-1.)
451            ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
452            pp=ppde.getSolution()
453            uu=loc(pp)
454            
455            # arguments for DcRes
456            current = 10.
457            sourceInfo = [ "sl0",  "sl1" ]
458            sampleTags = [ ("sr0", "sr1") ]
459    
460            sigmaPrimary=7.
461            phiPrimary=pp*current*sigma0/sigmaPrimary
462    
463            uuscale=1-current*sigma0/sigmaPrimary
464            delphi_in = [ (uu[1]-uu[0]) * uuscale]
465            
466            acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)
467    
468            self.assertTrue(Lsup(phiPrimary-acw.getPrimaryPotential()) < 1.e-10 * Lsup(acw.getPrimaryPotential()))
469    
470            SIGMA=10. # matches current        
471            args0=acw.getArguments(SIGMA)
472            p=args0[0]
473            u=args0[1]
474    
475            # true secondary potential
476            pps=pp-phiPrimary
477            self.assertTrue(Lsup(p-pps) < 1.e-6 * Lsup(pps))
478    
479    
480            # test return values at electrodes:
481            self.assertTrue(abs(u[0]-uu[0]*uuscale) < 1.e-6 * abs(uu[0]*uuscale))
482            self.assertTrue(abs(u[1]-uu[1]*uuscale) < 1.e-6 * abs(uu[1]*uuscale))
483    
484            # this sould be zero
485            dd=acw.getDefect(SIGMA, *args0)
486            self.assertTrue( dd >= 0.)
487            self.assertTrue( dd <= 1e-7 )
488    
489        def test_Differential2D(self):
490    
491            INC=0.001
492    
493            sigma0=1.
494            dx_tests=0.1
495            electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
496            from esys.finley import Rectangle
497            domain=Rectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
498            loc=Locator(domain,electrodes[2:])
499    
500            # arguments for DcRes
501            #current = 10.
502            sampleTags = [ ("sr0", "sr1") ]
503    
504            delphi_in = [ 0.05 ]
505    
506            sigmaPrimary=1
507            x=domain.getX()
508            phiPrimary=(x[0]-inf(x[0]))*(x[1]-inf(x[1]))*(x[0]-sup(x[0]))
509    
510            acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)
511    
512            #===========================================================================
513            x=Function(domain).getX()
514            SIGMA0=x[0]*x[1]+1
515            args0=acw.getArguments(SIGMA0)
516            d0=acw.getDefect(SIGMA0, *args0)
517            grad_d=acw.getGradient(SIGMA0, *args0)
518            
519            dS=exp(-(length(x-[0.5,0.5])/0.2)**2)
520            SIGMA1=SIGMA0+INC*dS
521            args1=acw.getArguments(SIGMA1)
522            d1=acw.getDefect(SIGMA1, *args1)
523            ref=abs((d1-d0)/INC)
524            self.assertTrue(abs((d1-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
525    
526            dS=-exp(-(length(x-[0.5,0.5])/0.2)**2)
527            SIGMA2=SIGMA0+INC*dS
528            args2=acw.getArguments(SIGMA2)
529            d2=acw.getDefect(SIGMA2, *args2)
530            ref=abs((d2-d0)/INC)
531            self.assertTrue(abs((d2-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
532    
533            dS=-1
534            SIGMA3=SIGMA0+INC*dS
535            args3=acw.getArguments(SIGMA3)
536            d3=acw.getDefect(SIGMA3, *args3)
537            ref=abs((d3-d0)/INC)
538            self.assertTrue(abs((d3-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
539    
540            dS=1
541            SIGMA4=SIGMA0+INC*dS
542            args4=acw.getArguments(SIGMA4)
543            d4=acw.getDefect(SIGMA4, *args4)
544            ref=abs((d4-d0)/INC)
545            self.assertTrue(abs((d4-d0)/INC-integrate(grad_d* dS)) < ref * 1.e-3)
546    
547  class TestIsostaticPressure(unittest.TestCase):  class TestIsostaticPressure(unittest.TestCase):
548      def test_all(self):      def test_all(self):

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

  ViewVC Help
Powered by ViewVC 1.1.26