/[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 1552 by gross, Thu May 8 08:52:41 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 2008 UTC
# Line 1  Line 1 
1  #  #
2  #   AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION  #   AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
3  #  #
4  #  #
5  #    step 1 rho*(v_star-v) = dt * (sigma'_ij,j-teta3*p,i+f_i)  #    step 1 rho*(v_star-v) = dt * (sigma'_ij,j-teta3*p,i+f_i)
6  #    step 2 dp=-dt*B*(v_j,j+teta1*v_star_j,j-dt*teta1*((1-teta3)*p_,jj+teta2*dp_,jj))  #    step 2 dp=-dt*B*(v_j,j+teta1*v_star_j,j-dt*teta1*((1-teta3)*p_,jj+teta2*dp_,jj))
7  #    step 3 rho*(v+-v) = -dt*((1-teta3)*p_,jj+teta2*dp_,jj)  #    step 3 rho*(v+-v) = -dt*((1-teta3)*p_,jj+teta2*dp_,jj)
8  #    step 4 p+=1/2(p+dp+abs(p+dp))  #    step 3b p+=1/2(p+dp+abs(p+dp))
9  #    step 4 sigma'i+_ij,j=f(v+,p+,...)  #    step 4 sigma'i+_ij,j=f(v+,p+,...)
10  #  #
11  #  #
12  from esys.escript import *  from esys.escript import *
13  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem  from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
14  from esys.finley import Rectangle  from esys.finley import Rectangle
15  iter     =   0  
16  nstep    =   3000  
17  w_step=max(int(nstep/50),1)*0+1  nel      =   20
18  mstep    =   5  H        =   0.5
19  H        =   0.5  L        =   1.0
20  eta      =   1.0       # shear viscosity  
21  ro       =   1.0  eta      =   1.0       # shear viscosity
22  g        =   1.00  ro       =   1.0
23  tauy0    =   0.9*ro*g*H/sqrt(3)       #0.11/(3*sqrt(3))  g        =   1.00
24  Pen      =   100  
25  vtop     =   0.01  alpha_w   =   1.00
26  dt       =   10**(-10)  alpha    =   1.00*1000000.
27  small    =   EPSILON  Pen=0.
28  alpha    =   1.00  B=100.
29  nel      =   20  
30  # len      =   sqrt(0.5/nel)  nstep    =   3000
31  toler    =   0.001  dt       =   1.
32  alphaw   =   1.00  small    =   EPSILON
33  L        =   1.0  w_step=max(int(nstep/50),1)*0+1
34  teta1    =    0.5  toler    =   0.001
35  teta2    =    0.5  teta1    =    0.5
36  teta3    =    0  # =0 split A; =1 split B  teta2    =    0.5
37  Etau=1000000000.  teta3    =    1  # =0 split A; =1 split B
38    
39  # create domain:  # create domain:
40  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)
41  x=dom.getX()  x=dom.getX()
42    
43    
44  momentumStep1=LinearPDESystem(dom) # A momentumStep1  momentumStep1=LinearPDESystem(dom)
45  momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])  momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) # fix x0=0 and x1=0
46  face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])  face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
47    
48                        # b   1 U1_2={0.0}  pressureStep2=LinearSinglePDE(dom)
49                        # b   1 U1_1={0,0}  pressureStep2.setReducedOrderOn()
50                                                   # b   4 U1_1={0.0}  pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
51    
52  pressureStep2=LinearSinglePDE(dom) # A [reduced] pressureStep2  
53  pressureStep2.setReducedOrderOn() #  
54  pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))  momentumStep3=LinearPDESystem(dom)
55                             # b  free U1={0.0} ?  momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])
56    #
57    #   initial values:
58    #
59  momentumStep3=LinearPDESystem(dom) # A momentumStep3  U=Vector(0.,Solution(dom))
60  momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3
61                        # b   1 U1_2={0.0}  # p=ro*g*(H-ReducedSolution(dom).getX()[1])
62                        # b   1 U1_1={0,0}  stress=Tensor(0.,Function(dom))
63                                                   # b   4 U1_1={0.0}  
64    t=dt
65    istep=0
66  # A deformation  while istep < nstep:
67  # e  V700= [grad] V100      istep=istep+1
68  # e  V901= V701+V704 +V101/(X1+small)      print "time step :",istep," t = ",t
69  # e  V601=sqrt(2*((V701)^2+(V704)^2+(V101/(X1+small))^2+\      r=Function(dom).getX()[0]
70  # (V702+V703)^2/2))      r_b=FunctionOnBoundary(dom).getX()[0]
71  # e  U1={abs(V901)/(V601+small)}      print " volume : ",integrate(r)
72        #
73        #  step 1:
74        #
75        # calculate normal
76  # A stress      n_d=dom.getNormal()
77  # e V700= [grad] V100      t_d=matrixmult(numarray.array([[0.,-1.],[1.,0]]),n_d)
78  # e  V901= V701+V704 +V101/(X1+small)      sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)
79  # e  V601=sqrt(2*((V701-V901)^2+(V704-V901)^2+(V101/(X1+small)-V901)^2+\      print " sigma_d =",inf(sigma_d),sup(sigma_d)
80  # (V702+V703)^2/2))  
81  # e  V602=((eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small)))      momentumStep1.setValue(D=ro*kronecker(dom),
82  # e V704=2*V602*V704-V401                             Y=ro*U+dt*((stress[:,0]-p*teta3*kronecker(dom)[:,0])/r+[0.,-ro*g]),
83  # b 3 [integrated] {V704}                             # Y=r*ro*U+dt*r*[0.,-ro*g],
84  #                             X=-dt*(stress-teta3*p*kronecker(dom)),
85  # A surface                             y=sigma_d*face_mask)
86  # b 3 [integrated] {1}      U_star=momentumStep1.getSolution()
87        #
88        #  step 2:
89  U=Vector(0.,Solution(dom)) # V100=0      #
90  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 # V401=ro*g*(L-X1)*(H-X2)/3      # U2=U+teta1*(U_star-U)
91  t=dt      U2=U+teta1*U_star
92  istep=0      gg2=grad(U2)
93  En=0      div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r
94    
95        grad_p=grad(p)
96  while istep < nstep:  
97      istep=istep+1      pressureStep2.setValue(A=r*dt*B*teta1*teta2/ro*dt*kronecker(dom),
98      print "time step :",istep," t = ",t                             D=r,                            
99                               Y=-dt*B*r*div_U2,
100      n_d=dom.getNormal()                             X=-r*B*dt**2/ro*teta1*(1-teta3)*grad_p)
101      t_d=matrixmult(numarray.array([[0.,1.],[-1.,0]]),n_d)      dp=pressureStep2.getSolution()
102        #
103      s=-En*inner(U,n_d)*face_mask      #  step 3:
104      nStressn=s*wherePositive(s)      #
105      u_boundary=Etau*inner(U,t_d)*face_mask      p2=(1-teta3)*p+teta2*dp
106      m=whereNegative(u_boundary-nStressn)      grad_p2=grad(p2)
107      nStressTau=nStressn*(1.-m)+u_boundary*m      momentumStep3.setValue(D=r*ro*kronecker(dom),
108      print nStressn                             Y=r*(ro*U_star-dt*teta2*grad_p2))
109      # print nStressn*n_d+nStressTau*t_d      U_new=momentumStep3.getSolution()
110      if istep == 20:      #
111          # print nStressTau      #   update:
112          saveVTK("stress.vtu",s=(nStressn*n_d+nStressTau*t_d))      #
113          1/0      p+=dp        
114      # print nStressTau      U=U_new
115        print " U:",inf(U),sup(U)
116      r=dom.getX()[0]      print " P:",inf(p),sup(p)
117            # A viscosity  
118      gg=grad(U) # V700= [grad] V100  
119      vol=gg[0,0]+gg[1,1]+U[0]/(r+small)    # V901= V701+V704 +V101/(X1+small)      p_pos=clip(p,small)
120      tau=sqrt(2*((gg[0,0]-vol/3)**2+(gg[1,1]-vol/3)**2+(U[0]/(r+small)-vol/3)**2+(gg[1,0]+gg[0,1])**2/2))      gg=grad(U)
121                # V601=sqrt(2*((V701-V901/3)^2+(V704-V901/3)^2+(V101/(X1+small)-V901/3)^2+(V702+V703)^2/2))      vol=gg[0,0]+gg[1,1]+U[0]/r  
122      m=whereNegative(eta*tau-alpha*p) # eta*V601<alpha*V401      gamma=sqrt(2*((gg[0,0]-vol/3)**2+(gg[1,1]-vol/3)**2+(U[0]/r-vol/3)**2+(gg[1,0]+gg[0,1])**2/2))
123      eta_d=m*eta+(1.-m)*alpha*p/(tau+small) #      m=whereNegative(eta*gamma-alpha*p_pos)
124                 # U1={(eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))}      eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)  
125            # viscosity 1 100 V100 200 V200 300 V300 400 V400      print " viscosity =",inf(eta_d),sup(eta_d)
126            # solve      stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
127      print " viscosity =",inf(eta_d),sup(eta_d) # show V801      #
128      stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))      # step size control:
129                 # solve momentum equationStep1      #
130                 # momentumStep1 1 100 V100 200 V200 300 V300 400 V400      len=inf(dom.getSize())
131      momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i      dt1=inf(dom.getSize()/(length(U)+small))
132                             Y=stress[:,0]/(r+small)+[0.,-ro*g],           # -{0.0,ro*g}_i      dt2=inf(0.5*ro*(len**2)/eta_d)
133                             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}      dt=dt1*dt2/(dt1+dt2)
134                             y=-(nStressn*n_d+nStressTau*t_d))      print " new step size = ",dt
135      dU2=momentumStep1.getSolution()      #
136            # solve      #  update geometry
137      #                     dU2-> V100      #
138      #                     U -> V200      dom.setX(dom.getX()+U*dt)
139      #                     p -> V500      t=t+dt
140        if (istep-1)%w_step==0:saveVTK("u.%d.xml"%((istep-1)/w_step),p=p,eta=eta_d,U=U_star,U_star=U_star,gamma=gamma)
141      # pressureStep2 1 100 V100 200 V200 300 V300 400 V500          if istep == 3: 1/0
       
     gg2=grad(dU2)  
     div_dU2=gg2[0,0]+gg2[1,1]+dU2[0]/(r+small)  
     pressureStep2.setValue(A=r*teta1*teta2*dt**2*kronecker(dom),  # -D_j{teta1*teta2*dt^2}D_jU1  
                            D=r*ro,                            # {ro/Pen}U1  
                            Y=-r*(ro*dt*vol+teta1*ro*dt*div_dU2))  # -{ro*dt}D_j{V200}_j-{teta1*ro*dt}D_j{V100}_j  
     # solve  
     dp=pressureStep2.getSolution()  
                    # dp -> V100  
                    # dU2 -> V200  
                    # U  -> V300  
                    # p -> V600  
     p+=dp                  # V601=V601+V101  
                            # V701=V101 : dp -> V700  
     p=(p+abs(p))/2+small # V601=cut V601  
                          # V601=V601+small  
     print " pressure increment :",inf(dp),sup(dp) # show V601  
     print " pressure range :",inf(p),sup(p) # show V601  
                 # popp  
                 #   dU2 -> V100  
                 #   U  -> V200  
                 #   p -> V500  
                 #   dp -> V600  
   
     # momentumStep3 1 100 V100 200 V200 300 V300 400 V600  
     A2=Tensor4(0.0,Function(dom))  
     for i in range(dom.getDim()):  
        for j in range(dom.getDim()):  
           A2[i,i,j,j]+=1  
     momentumStep3.setValue(A=Pen*A2,  # -D_i{Pen}D_jU1_j  
                            D=ro/dt*kronecker(dom), # {ro/dt}U1_i  
                            Y=(ro/dt*(dU2+U)-teta2*grad(dp)))  
                            # Y=r*(ro/dt*(dU2+U)-teta2*dp*[1.,0]),       # {ro/dt}{V200}_i+{ro/dt}{V100}_i  
                            # X=r*teta2*dp**kronecker(dom))            # -D_i{teta2*V401}  
     U, U_old=momentumStep3.getSolution(), U  
                 #   U -> V100  
                 #   dU2 -> V200  
                 #   U_old  -> V300  
                 #   p -> V600  
                 #   dp -> V700  
            # V400=V300  
            # V300=V100  
            # popp  
            # popp  
                 #   U  -> V100  
                 #   U_old  -> V200  
                 #   p -> V400  
                 #   dp -> V500  
            # deformation  
            # solve  
            # show ratio  
            # show V101  
            # popp  
   
     print " new U:",inf(U),sup(U) # show v100  
     print " old U",inf(U_old),sup(U_old) # show v200  
     print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200  
                                                     # test0=L1 V100  
                                                     # show test/test0  
                      # black  
                      # arrow  
         
     vmax=Lsup(U)  # vmax1=abs(max V100)  
                            # vmax2=abs(min V100)  
                            # show vmax1  
                            # show vmax2  
         
                            # vmax=vmax1  
                            # if vmax < vmax2  
                            # vmax=vmax2  
                            # endif  
     print " max velocity = ",vmax # show vmax  
     len=inf(dom.getSize())  
     dt1=len/vmax  # Courant condition  
     dt2=0.5*ro*(len**2)/eta  
     dt=dt1*dt2/(dt1+dt2)  
     En=1/dt  
     t=t+dt  
     print " time , time step ",t,dt # show t  
                                         # show dt  
           
   
     dom.setX(dom.getX()+U*dt) # mapm 500  
                             # V500=V500+V100*dt  
                             # V600=V500  
                             # mapm 500  
   
     print " volume : ",integrate(r)# vol=Out  
                                           # show vol  
                 #   dU  -> V100  
                 #   U  -> V200  
                 #   p -> V400  
                 #   dp -> V600  
                # clear  
                # prim  
     # saveVTK("u.%d.xml"%(istep),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp)  
     if (istep-1)%w_step==0:saveVTK("u.%d.xml"%((istep-1)/w_step),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp)  

Legend:
Removed from v.1552  
changed lines
  Added in v.1562

  ViewVC Help
Powered by ViewVC 1.1.26