/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

Diff of /trunk/finley/test/python/run_models.py

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

revision 2300 by gross, Wed Feb 11 06:48:28 2009 UTC revision 2301 by gross, Thu Mar 12 01:35:26 2009 UTC
# Line 37  VERBOSE=False  # or True Line 37  VERBOSE=False  # or True
37  DETAIL_VERBOSE=False  DETAIL_VERBOSE=False
38    
39  from esys.escript import *  from esys.escript import *
40  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian, PowerLaw
41  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
42    
43  from esys.escript.models import Mountains  from esys.escript.models import Mountains
# Line 742  class Test_Mountains2D(unittest.TestCase Line 742  class Test_Mountains2D(unittest.TestCase
742         error_int=integrate(Z)         error_int=integrate(Z)
743         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")         self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
744                
745    
746    
747    class Test_Rheologies(unittest.TestCase):
748    
749         """
750         this is the program used to generate the powerlaw tests:
751    
752         TAU_Y=100.
753         N=10
754         M=5
755    
756         def getE(tau):
757             if tau<=TAU_Y:
758               return 1./(0.5+20*sqrt(tau))
759             else:
760               raise ValueError,"out of range."
761         tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
762         e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
763         e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
764    
765         print tau
766         print e
767         """
768         TOL=1.e-8
769         def test_PowerLaw_Init(self):
770             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
771    
772             self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
773             self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
774             self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
775             self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
776             self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
777    
778             self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
779             self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
780             pl.setDruckerPragerLaw(tau_Y=10,friction=3)
781             self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
782             self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
783    
784             self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
785             pl.setElasticShearModulus(1000)
786             self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
787    
788             e=pl.getEtaN()
789             self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
790             self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
791             self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
792             self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
793             self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
794             self.failUnlessRaises(ValueError, pl.getEtaN, 3)
795    
796             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
797             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
798             self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
799    
800             pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
801             self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
802             self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
803             self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
804    
805             pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
806             self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
807             self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
808             self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
809    
810             self.failUnlessRaises(ValueError,pl.getPower,-1)
811             self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
812             self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
813             self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
814             self.failUnlessRaises(ValueError,pl.getPower,3)
815    
816             self.failUnlessRaises(ValueError,pl.getTauT,-1)
817             self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
818             self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
819             self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
820             self.failUnlessRaises(ValueError,pl.getTauT,3)
821    
822             self.failUnlessRaises(ValueError,pl.getEtaN,-1)
823             self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
824             self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
825             self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
826             self.failUnlessRaises(ValueError,pl.getEtaN,3)
827    
828         def checkResult(self,id,gamma_dot_, eta, tau_ref):
829             self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
830             error=abs(gamma_dot_*eta-tau_ref)
831             self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id))
832            
833         def test_PowerLaw_Linear(self):
834             taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
835             gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0]
836             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
837             pl.setDruckerPragerLaw(tau_Y=100.)
838             pl.setPowerLaw(eta_N=2.)
839             pl.setEtaTolerance(self.TOL)
840             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
841            
842         def test_PowerLaw_QuadLarge(self):
843             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
844             gamma_dot_s=[0.0, 637.45553203367592, 1798.8543819998317, 3301.3353450309969, 5079.6442562694074, 7096.067811865476, 9325.1600308977995, 11748.240371477059, 14350.835055998656, 17121.29936490925, 20050.0, 22055.0, 24060.0, 26065.0, 28070.0, 30075.0]
845             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
846             pl.setDruckerPragerLaw(tau_Y=100.)
847             pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
848             pl.setEtaTolerance(self.TOL)
849             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
850    
851         def test_PowerLaw_QuadSmall(self):
852             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
853             gamma_dot_s=[0.0, 5.632455532033676, 11.788854381999831, 18.286335345030995, 25.059644256269408, 32.071067811865476, 39.295160030897804, 46.713240371477056, 54.310835055998652, 62.076299364909254, 70.0, 77.0, 84.0, 91.0, 98.0, 105.0]
854             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
855             pl.setDruckerPragerLaw(tau_Y=100.)
856             pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
857             pl.setEtaTolerance(self.TOL)
858             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
859    
860         def test_PowerLaw_CubeLarge(self):
861             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
862             gamma_dot_s=[0.0, 51.415888336127786, 157.36125994561547, 304.6468153816889, 487.84283811405851, 703.60440414872664, 949.57131887226353, 1223.9494765692673, 1525.3084267560891, 1852.4689652218574, 2204.4346900318833, 2424.8781590350718, 2645.3216280382603, 2865.7650970414484, 3086.2085660446369, 3306.6520350478249]
863             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
864             pl.setDruckerPragerLaw(tau_Y=100.)
865             pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
866             pl.setEtaTolerance(self.TOL)
867             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
868    
869         def test_PowerLaw_CubeSmall(self):
870             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
871             gamma_dot_s=[0.0, 5.4641588833612778, 11.473612599456157, 17.89646815381689, 24.678428381140588, 31.786044041487269, 39.195713188722635, 46.889494765692675, 54.853084267560895, 63.074689652218574, 71.544346900318828, 78.698781590350706, 85.853216280382583, 93.007650970414474, 100.16208566044635, 107.316520350478]
872             pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
873             pl.setDruckerPragerLaw(tau_Y=100.)
874             pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
875             pl.setEtaTolerance(self.TOL)
876             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
877    
878         def test_PowerLaw_QuadLarge_CubeLarge(self):
879             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
880             gamma_dot_s=[0.0, 683.87142036980367, 1946.2156419454475, 3590.982160412686, 5547.4870943834667, 7774.6722160142008, 10244.731349770063, 12937.189848046326, 15836.143482754744, 18928.768330131104, 22204.434690031885, 24424.878159035074, 26645.321628038262, 28865.765097041451, 31086.208566044639, 33306.652035047824]
881             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
882             pl.setDruckerPragerLaw(tau_Y=100.)
883             pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
884             pl.setEtaTolerance(self.TOL)
885             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
886    
887         def test_PowerLaw_QuadLarge_CubeSmall(self):
888             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
889             gamma_dot_s=[0.0, 637.9196909170372, 1800.3279945992881, 3304.2318131848137, 5084.3226846505486, 7102.853855906963, 9334.3557440865225, 11760.129866242751, 14365.688140266215, 17139.374054561471, 20071.544346900318, 22078.698781590349, 24085.853216280382, 26093.007650970416, 28100.162085660446, 30107.316520350476]
890             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
891             pl.setDruckerPragerLaw(tau_Y=100.)
892             pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
893             pl.setEtaTolerance(self.TOL)
894             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
895    
896         def test_PowerLaw_QuadSmall_CubeLarge(self):
897             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
898             gamma_dot_s=[0.0, 52.04834386816146, 159.15011432761528, 307.93315072671987, 492.9024823703279, 710.67547196059206, 958.86647890316135, 1235.6627169407443, 1539.6192618120876, 1869.5452645867665, 2224.4346900318833, 2446.8781590350718, 2669.3216280382603, 2891.7650970414484, 3114.2085660446369, 3336.6520350478249]
899    
900             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
901             pl.setDruckerPragerLaw(tau_Y=100.)
902             pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
903             pl.setEtaTolerance(self.TOL)
904             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
905    
906         def test_PowerLaw_QuadSmall_CubeSmall(self):
907             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
908             gamma_dot_s=[0.0, 6.0966144153949529, 13.262466981455987, 21.182803498847885, 29.738072637409996, 38.857111853352741, 48.49087321962044, 58.602735137169738, 69.16391932355954, 80.150989017127827, 91.544346900318828, 100.69878159035071, 109.85321628038258, 119.00765097041447, 128.16208566044637, 137.31652035047824]
909             pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
910             pl.setDruckerPragerLaw(tau_Y=100.)
911             pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
912             pl.setEtaTolerance(self.TOL)
913             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
914    
915         def test_PowerLaw_withShear(self):
916             taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
917             gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
918             pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
919             pl.setDruckerPragerLaw(tau_Y=100.)
920             pl.setPowerLaw(eta_N=2.)
921             pl.setElasticShearModulus(3.)
922             dt=1./3.
923             pl.setEtaTolerance(self.TOL)
924             self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
925             for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
926    
927  if __name__ == '__main__':  if __name__ == '__main__':
928     suite = unittest.TestSuite()     suite = unittest.TestSuite()
929        
# Line 751  if __name__ == '__main__': Line 933  if __name__ == '__main__':
933     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
934     suite.addTest(unittest.makeSuite(Test_Mountains3D))     suite.addTest(unittest.makeSuite(Test_Mountains3D))
935     suite.addTest(unittest.makeSuite(Test_Mountains2D))     suite.addTest(unittest.makeSuite(Test_Mountains2D))
936       suite.addTest(unittest.makeSuite(Test_Rheologies))
937     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
938     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
939    
Legend:
Removed from v.2300  
changed lines
  Added in v.2301
  ViewVC Help
Powered by ViewVC 1.1.26