/[escript]/trunk/finley/test/python/axisymm-splitB.py
ViewVC logotype

Diff of /trunk/finley/test/python/axisymm-splitB.py

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

revision 1392 by gross, Mon Jan 21 06:20:54 2008 UTC revision 1417 by gross, Mon Feb 25 04:45:48 2008 UTC
# Line 6  from esys.escript.linearPDEs import Line Line 6  from esys.escript.linearPDEs import Line
6  from esys.finley import Rectangle  from esys.finley import Rectangle
7  iter     =   0  iter     =   0
8  nstep    =   3000  nstep    =   3000
9    w_step=max(int(nstep/50),1)*0+1
10  mstep    =   5  mstep    =   5
11  H        =   0.5  H        =   0.5
12  eta      =   1.0       # shear viscosity  eta      =   1.0       # shear viscosity
# Line 24  alphaw   =   1.00 Line 25  alphaw   =   1.00
25  L        =   1.0  L        =   1.0
26  teta1    =    0.5  teta1    =    0.5
27  teta2    =    0.5  teta2    =    0.5
28    Etau=1000000000.
29    
30  dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)  dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)
31  x=dom.getX()  x=dom.getX()
32    
33  momentumStep1=LinearPDESystem(dom) # A momentumStep1  momentumStep1=LinearPDESystem(dom) # A momentumStep1
34  momentumStep1.setValue(q=whereZero(x[1])*[1.,1.]+whereZero(x[0])*[1.,0.])  momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
35    face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
36    
37                        # b   1 U1_2={0.0}                        # b   1 U1_2={0.0}
38                        # b   1 U1_1={0,0}                        # b   1 U1_1={0,0}
39                                                   # b   4 U1_1={0.0}                                                   # b   4 U1_1={0.0}
# Line 42  pressureStep2.setValue(q=whereZero(x[0]- Line 46  pressureStep2.setValue(q=whereZero(x[0]-
46    
47    
48  momentumStep3=LinearPDESystem(dom) # A momentumStep3  momentumStep3=LinearPDESystem(dom) # A momentumStep3
49  momentumStep3.setValue(q=whereZero(x[1])*[1.,1.]+whereZero(x[0])*[1.,0.])  momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
50                        # b   1 U1_2={0.0}                        # b   1 U1_2={0.0}
51                        # b   1 U1_1={0,0}                        # b   1 U1_1={0,0}
52                                                   # b   4 U1_1={0.0}                                                   # b   4 U1_1={0.0}
# Line 75  U=Vector(0.,Solution(dom)) # V100=0 Line 79  U=Vector(0.,Solution(dom)) # V100=0
79  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 # V401=ro*g*(L-X1)*(H-X2)/3  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 # V401=ro*g*(L-X1)*(H-X2)/3
80  t=dt  t=dt
81  istep=0  istep=0
82  w_step=max(int(nstep/50),1)  En=0
83    
84    
85  while istep < nstep:  while istep < nstep:
86      istep=istep+1      istep=istep+1
87      print "time step :",istep," t = ",t      print "time step :",istep," t = ",t
88    
89        n_d=dom.getNormal()
90        t_d=matrixmult(numarray.array([[0.,1.],[-1.,0]]),n_d)
91    
92        s=-En*inner(U,n_d)*face_mask
93        nStressn=s*wherePositive(s)
94        u_boundary=Etau*inner(U,t_d)*face_mask
95        m=whereNegative(u_boundary-nStressn)
96        nStressTau=nStressn*(1.-m)+u_boundary*m
97        print nStressn
98        # print nStressn*n_d+nStressTau*t_d
99        if istep == 20:
100            # print nStressTau
101            saveVTK("stress.vtu",s=(nStressn*n_d+nStressTau*t_d))
102            1/0
103        # print nStressTau
104    
105      r=dom.getX()[0]      r=dom.getX()[0]
106            # A viscosity            # A viscosity
107      gg=grad(U) # V700= [grad] V100      gg=grad(U) # V700= [grad] V100
# Line 98  while istep < nstep: Line 119  while istep < nstep:
119                 # momentumStep1 1 100 V100 200 V200 300 V300 400 V400                 # momentumStep1 1 100 V100 200 V200 300 V300 400 V400
120      momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i      momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i
121                             Y=stress[:,0]/(r+small)+[0.,-ro*g],           # -{0.0,ro*g}_i                             Y=stress[:,0]/(r+small)+[0.,-ro*g],           # -{0.0,ro*g}_i
122                             X=-(stress+p*kronecker(dom))) # -D_j{V602}D_j{V100}_i-D_i{V602}D_j{V100}_j+D_i{(2/3)*V602}D_j{V100}_j-D_i{V401}                             X=-(stress+p*kronecker(dom)), # -D_j{V602}D_j{V100}_i-D_i{V602}D_j{V100}_j+D_i{(2/3)*V602}D_j{V100}_j-D_i{V401}
123                               y=-(nStressn*n_d+nStressTau*t_d))
124      dU2=momentumStep1.getSolution()      dU2=momentumStep1.getSolution()
125            # solve            # solve
126      #                     dU2-> V100      #                     dU2-> V100
# Line 159  while istep < nstep: Line 181  while istep < nstep:
181             # show ratio             # show ratio
182             # show V101             # show V101
183             # popp             # popp
184    
185      print " new U:",inf(U),sup(U) # show v100      print " new U:",inf(U),sup(U) # show v100
186      print " old U",inf(U_old),sup(U_old) # show v200      print " old U",inf(U_old),sup(U_old) # show v200
187      print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200      print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200
# Line 181  while istep < nstep: Line 204  while istep < nstep:
204      dt1=len/vmax  # Courant condition      dt1=len/vmax  # Courant condition
205      dt2=0.5*ro*(len**2)/eta      dt2=0.5*ro*(len**2)/eta
206      dt=dt1*dt2/(dt1+dt2)      dt=dt1*dt2/(dt1+dt2)
207        En=1/dt
208      t=t+dt      t=t+dt
209      print " time , time step ",t,dt # show t      print " time , time step ",t,dt # show t
210                                          # show dt                                          # show dt

Legend:
Removed from v.1392  
changed lines
  Added in v.1417

  ViewVC Help
Powered by ViewVC 1.1.26