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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1392 - (show annotations)
Mon Jan 21 06:20:54 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 7576 byte(s)
axissymmetric collaps added.
1 #
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