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

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
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)),
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
93  En=0      div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r
94
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,
101      t_d=matrixmult(numarray.array([[0.,1.],[-1.,0]]),n_d)      dp=pressureStep2.getSolution()
102        #
104      nStressn=s*wherePositive(s)      #
107      nStressTau=nStressn*(1.-m)+u_boundary*m      momentumStep3.setValue(D=r*ro*kronecker(dom),
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
119      vol=gg[0,0]+gg[1,1]+U[0]/(r+small)    # V901= V701+V704 +V101/(X1+small)      p_pos=clip(p,small)
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

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=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