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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1392 - (hide annotations)
Mon Jan 21 06:20:54 2008 UTC (13 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 7576 byte(s)
axissymmetric collaps added.
1 gross 1392 #
2     # AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
3     #
4     from esys.escript import *
5     from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
6     from esys.finley import Rectangle
7     iter = 0
8     nstep = 3000
9     mstep = 5
10     H = 0.5
11     eta = 1.0 # shear viscosity
12     ro = 1.0
13     g = 1.00
14     tauy0 = 0.9*ro*g*H/sqrt(3) #0.11/(3*sqrt(3))
15     Pen = 100
16     vtop = 0.01
17     dt = 10**(-10)
18     small = EPSILON
19     alpha = 1.00
20     nel = 20
21     # len = sqrt(0.5/nel)
22     toler = 0.001
23     alphaw = 1.00
24     L = 1.0
25     teta1 = 0.5
26     teta2 = 0.5
27    
28     dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)
29     x=dom.getX()
30    
31     momentumStep1=LinearPDESystem(dom) # A momentumStep1
32     momentumStep1.setValue(q=whereZero(x[1])*[1.,1.]+whereZero(x[0])*[1.,0.])
33     # b 1 U1_2={0.0}
34     # b 1 U1_1={0,0}
35     # b 4 U1_1={0.0}
36    
37     pressureStep2=LinearSinglePDE(dom) # A [reduced] pressureStep2
38     pressureStep2.setReducedOrderOn() #
39     pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
40     # b free U1={0.0} ?
41    
42    
43    
44     momentumStep3=LinearPDESystem(dom) # A momentumStep3
45     momentumStep3.setValue(q=whereZero(x[1])*[1.,1.]+whereZero(x[0])*[1.,0.])
46     # b 1 U1_2={0.0}
47     # b 1 U1_1={0,0}
48     # b 4 U1_1={0.0}
49    
50    
51     # A deformation
52     # e V700= [grad] V100
53     # e V901= V701+V704 +V101/(X1+small)
54     # e V601=sqrt(2*((V701)^2+(V704)^2+(V101/(X1+small))^2+\
55     # (V702+V703)^2/2))
56     # e U1={abs(V901)/(V601+small)}
57    
58    
59    
60    
61     # A stress
62     # e V700= [grad] V100
63     # e V901= V701+V704 +V101/(X1+small)
64     # e V601=sqrt(2*((V701-V901)^2+(V704-V901)^2+(V101/(X1+small)-V901)^2+\
65     # (V702+V703)^2/2))
66     # e V602=((eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small)))
67     # e V704=2*V602*V704-V401
68     # b 3 [integrated] {V704}
69     #
70     # A surface
71     # b 3 [integrated] {1}
72    
73    
74     U=Vector(0.,Solution(dom)) # V100=0
75     p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 # V401=ro*g*(L-X1)*(H-X2)/3
76     t=dt
77     istep=0
78     w_step=max(int(nstep/50),1)
79    
80    
81     while istep < nstep:
82     istep=istep+1
83     print "time step :",istep," t = ",t
84     r=dom.getX()[0]
85     # A viscosity
86     gg=grad(U) # V700= [grad] V100
87     vol=gg[0,0]+gg[1,1]+U[0]/(r+small) # V901= V701+V704 +V101/(X1+small)
88     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))
89     # V601=sqrt(2*((V701-V901/3)^2+(V704-V901/3)^2+(V101/(X1+small)-V901/3)^2+(V702+V703)^2/2))
90     m=whereNegative(eta*tau-alpha*p) # eta*V601<alpha*V401
91     eta_d=m*eta+(1.-m)*alpha*p/(tau+small) #
92     # U1={(eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))}
93     # viscosity 1 100 V100 200 V200 300 V300 400 V400
94     # solve
95     print " viscosity =",inf(eta_d),sup(eta_d) # show V801
96     stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
97     # solve momentum equationStep1
98     # momentumStep1 1 100 V100 200 V200 300 V300 400 V400
99     momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i
100     Y=stress[:,0]/(r+small)+[0.,-ro*g], # -{0.0,ro*g}_i
101     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}
102     dU2=momentumStep1.getSolution()
103     # solve
104     # dU2-> V100
105     # U -> V200
106     # p -> V500
107    
108     # pressureStep2 1 100 V100 200 V200 300 V300 400 V500
109    
110     gg2=grad(dU2)
111     div_dU2=gg2[0,0]+gg2[1,1]+dU2[0]/(r+small)
112     pressureStep2.setValue(A=r*teta1*teta2*dt**2*kronecker(dom), # -D_j{teta1*teta2*dt^2}D_jU1
113     D=r*ro, # {ro/Pen}U1
114     Y=-r*(ro*dt*vol+teta1*ro*dt*div_dU2)) # -{ro*dt}D_j{V200}_j-{teta1*ro*dt}D_j{V100}_j
115     # solve
116     dp=pressureStep2.getSolution()
117     # dp -> V100
118     # dU2 -> V200
119     # U -> V300
120     # p -> V600
121     p+=dp # V601=V601+V101
122     # V701=V101 : dp -> V700
123     p=(p+abs(p))/2+small # V601=cut V601
124     # V601=V601+small
125     print " pressure increment :",inf(dp),sup(dp) # show V601
126     print " pressure range :",inf(p),sup(p) # show V601
127     # popp
128     # dU2 -> V100
129     # U -> V200
130     # p -> V500
131     # dp -> V600
132    
133     # momentumStep3 1 100 V100 200 V200 300 V300 400 V600
134     A2=Tensor4(0.0,Function(dom))
135     for i in range(dom.getDim()):
136     for j in range(dom.getDim()):
137     A2[i,i,j,j]+=1
138     momentumStep3.setValue(A=Pen*A2, # -D_i{Pen}D_jU1_j
139     D=ro/dt*kronecker(dom), # {ro/dt}U1_i
140     Y=(ro/dt*(dU2+U)-teta2*grad(dp)))
141     # Y=r*(ro/dt*(dU2+U)-teta2*dp*[1.,0]), # {ro/dt}{V200}_i+{ro/dt}{V100}_i
142     # X=r*teta2*dp**kronecker(dom)) # -D_i{teta2*V401}
143     U, U_old=momentumStep3.getSolution(), U
144     # U -> V100
145     # dU2 -> V200
146     # U_old -> V300
147     # p -> V600
148     # dp -> V700
149     # V400=V300
150     # V300=V100
151     # popp
152     # popp
153     # U -> V100
154     # U_old -> V200
155     # p -> V400
156     # dp -> V500
157     # deformation
158     # solve
159     # show ratio
160     # show V101
161     # popp
162     print " new U:",inf(U),sup(U) # show v100
163     print " old U",inf(U_old),sup(U_old) # show v200
164     print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200
165     # test0=L1 V100
166     # show test/test0
167     # black
168     # arrow
169    
170     vmax=Lsup(U) # vmax1=abs(max V100)
171     # vmax2=abs(min V100)
172     # show vmax1
173     # show vmax2
174    
175     # vmax=vmax1
176     # if vmax < vmax2
177     # vmax=vmax2
178     # endif
179     print " max velocity = ",vmax # show vmax
180     len=inf(dom.getSize())
181     dt1=len/vmax # Courant condition
182     dt2=0.5*ro*(len**2)/eta
183     dt=dt1*dt2/(dt1+dt2)
184     t=t+dt
185     print " time , time step ",t,dt # show t
186     # show dt
187    
188    
189     dom.setX(dom.getX()+U*dt) # mapm 500
190     # V500=V500+V100*dt
191     # V600=V500
192     # mapm 500
193    
194     print " volume : ",integrate(r)# vol=Out
195     # show vol
196     # dU -> V100
197     # U -> V200
198     # p -> V400
199     # dp -> V600
200     # clear
201     # prim
202     # saveVTK("u.%d.xml"%(istep),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp)
203     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)

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26