/[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 2431 by gross, Wed Apr 15 06:46:01 2009 UTC revision 2432 by gross, Wed May 20 06:06:20 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, PowerLaw  from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian
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 924  class Test_Rheologies(unittest.TestCase) Line 924  class Test_Rheologies(unittest.TestCase)
924           self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])           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])           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    class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
928       TOL=1.e-6
929       VERBOSE=False
930       A=1.
931       P_max=100
932       NE=2*getMPISizeWorld()
933       tau_Y=10.
934       N_dt=10
935    
936       # material parameter:
937       tau_1=5.
938       tau_2=5.
939       eta_0=100.
940       eta_1=50.
941       eta_2=400.
942       N_1=2.
943       N_2=3.
944       def getReference(self, t):
945    
946          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
947          x=self.dom.getX()
948    
949          s_00=min(self.A*t,B)
950          tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
951          inv_eta= 1./self.eta_0 + 1./self.eta_1*(tau/self.tau_1)**(1.-1./self.N_1) + 1./self.eta_2*(tau/self.tau_2)**(1.-1./self.N_2)
952    
953          alpha=0.5*inv_eta*s_00
954          if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
955          u_ref=x*alpha
956          u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
957          sigma_ref=kronecker(self.dom)*s_00
958          sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
959    
960          p_ref=self.P_max
961          for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
962          p_ref-=integrate(p_ref)/vol(self.dom)
963          return u_ref, sigma_ref, p_ref
964    
965       def runIt(self, free=None):
966          x=self.dom.getX()
967          B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
968          dt=B/int(self.N_dt/2)
969          if self.VERBOSE: print "dt =",dt
970          if self.latestart:
971              t=dt
972          else:
973              t=0
974          v,s,p=self.getReference(t)
975    
976          mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
977          mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
978          mod.setElasticShearModulus(self.mu)
979          mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2],  [1.,self.N_1,self.N_2])
980          mod.setTolerance(self.TOL)
981          mod.setFlowSubTolerance(self.TOL**2)
982          mod.setFlowTolerance(self.TOL)
983    
984          BF=Vector(self.P_max,Function(self.dom))
985          for d in range(self.dom.getDim()):
986              for d2 in range(self.dom.getDim()):
987                  if d!=d2: BF[d]*=x[d2]
988          v_mask=Vector(0,Solution(self.dom))
989          if free==None:
990             for d in range(self.dom.getDim()):
991                v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
992          else:
993             for d in range(self.dom.getDim()):
994                if d == self.dom.getDim()-1:
995                   v_mask[d]=whereZero(x[d]-1.)
996                else:
997                   v_mask[d]=whereZero(x[d])
998          mod.setExternals(F=BF,fixed_v_mask=v_mask)
999          
1000          n=self.dom.getNormal()
1001          N_t=0
1002          errors=[]
1003          while N_t < self.N_dt:
1004             t_ref=t+dt
1005             v_ref, s_ref,p_ref=self.getReference(t_ref)
1006             mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
1007             mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
1008             self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
1009             t+=dt
1010             N_t+=1
1011    
1012       def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1013             p=mod.getPressure()
1014             p-=integrate(p)/vol(self.dom)
1015             error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1016             error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1017             error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1018             error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1019             if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1020             self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1021             self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1022             self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1023             self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1024       def tearDown(self):
1025            del self.dom
1026    
1027       def test_D2_Fixed_MuNone_LateStart(self):
1028           self.dom = Rectangle(self.NE,self.NE,order=2)
1029           self.mu=None
1030           self.latestart=True
1031           self.runIt()
1032       def test_D2_Fixed_Mu_LateStart(self):
1033           self.dom = Rectangle(self.NE,self.NE,order=2)
1034           self.mu=555.
1035           self.latestart=True
1036           self.runIt()
1037       def test_D2_Fixed_MuNone(self):
1038           self.dom = Rectangle(self.NE,self.NE,order=2)
1039           self.mu=None
1040           self.latestart=False
1041           self.runIt()
1042       def test_D2_Fixed_Mu(self):
1043           self.dom = Rectangle(self.NE,self.NE,order=2)
1044           self.mu=555.
1045           self.latestart=False
1046           self.runIt()
1047       def test_D2_Free_MuNone_LateStart(self):
1048           self.dom = Rectangle(self.NE,self.NE,order=2)
1049           self.mu=None
1050           self.latestart=True
1051           self.runIt(free=0)
1052       def test_D2_Free_Mu_LateStart(self):
1053           self.dom = Rectangle(self.NE,self.NE,order=2)
1054           self.mu=555.
1055           self.latestart=True
1056           self.runIt(free=0)
1057       def test_D2_Free_MuNone(self):
1058           self.dom = Rectangle(self.NE,self.NE,order=2)
1059           self.mu=None
1060           self.latestart=False
1061           self.runIt(free=0)
1062       def test_D2_Free_Mu(self):
1063           self.dom = Rectangle(self.NE,self.NE,order=2)
1064           self.mu=555.
1065           self.latestart=False
1066           self.runIt(free=0)
1067    
1068       def test_D3_Fixed_MuNone_LateStart(self):
1069           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1070           self.mu=None
1071           self.latestart=True
1072           self.runIt()
1073       def test_D3_Fixed_Mu_LateStart(self):
1074           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1075           self.mu=555.
1076           self.latestart=True
1077           self.runIt()
1078       def test_D3_Fixed_MuNone(self):
1079           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1080           self.mu=None
1081           self.latestart=False
1082           self.runIt()
1083       def test_D3_Fixed_Mu(self):
1084           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1085           self.mu=555.
1086           self.latestart=False
1087           self.runIt()
1088       def test_D3_Free_MuNone_LateStart(self):
1089           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1090           self.mu=None
1091           self.latestart=True
1092           self.runIt(free=0)
1093       def test_D3_Free_Mu_LateStart(self):
1094           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1095           self.mu=555.
1096           self.latestart=True
1097           self.runIt(free=0)
1098       def test_D3_Free_MuNone(self):
1099           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1100           self.mu=None
1101           self.latestart=False
1102           self.runIt(free=0)
1103       def test_D3_Free_Mu(self):
1104           self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1105           self.mu=555.
1106           self.latestart=False
1107           self.runIt(free=0)
1108    
1109  if __name__ == '__main__':  if __name__ == '__main__':
1110     suite = unittest.TestSuite()     suite = unittest.TestSuite()
     
1111     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1112     suite.addTest(unittest.makeSuite(Test_Darcy3D))     suite.addTest(unittest.makeSuite(Test_Darcy3D))
1113     suite.addTest(unittest.makeSuite(Test_Darcy2D))     suite.addTest(unittest.makeSuite(Test_Darcy2D))
# Line 934  if __name__ == '__main__': Line 1115  if __name__ == '__main__':
1115     suite.addTest(unittest.makeSuite(Test_Mountains3D))     suite.addTest(unittest.makeSuite(Test_Mountains3D))
1116     suite.addTest(unittest.makeSuite(Test_Mountains2D))     suite.addTest(unittest.makeSuite(Test_Mountains2D))
1117     suite.addTest(unittest.makeSuite(Test_Rheologies))     suite.addTest(unittest.makeSuite(Test_Rheologies))
1118       suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1119     s=unittest.TextTestRunner(verbosity=2).run(suite)     s=unittest.TextTestRunner(verbosity=2).run(suite)
1120     if not s.wasSuccessful(): sys.exit(1)     if not s.wasSuccessful(): sys.exit(1)
1121    

Legend:
Removed from v.2431  
changed lines
  Added in v.2432

  ViewVC Help
Powered by ViewVC 1.1.26