/[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 1417 - (show annotations)
Mon Feb 25 04:45:48 2008 UTC (13 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 8245 byte(s)
some more work on the transport solver.
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 w_step=max(int(nstep/50),1)*0+1
10 mstep = 5
11 H = 0.5
12 eta = 1.0 # shear viscosity
13 ro = 1.0
14 g = 1.00
15 tauy0 = 0.9*ro*g*H/sqrt(3) #0.11/(3*sqrt(3))
16 Pen = 100
17 vtop = 0.01
18 dt = 10**(-10)
19 small = EPSILON
20 alpha = 1.00
21 nel = 20
22 # len = sqrt(0.5/nel)
23 toler = 0.001
24 alphaw = 1.00
25 L = 1.0
26 teta1 = 0.5
27 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)
31 x=dom.getX()
32
33 momentumStep1=LinearPDESystem(dom) # A momentumStep1
34 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}
38 # b 1 U1_1={0,0}
39 # b 4 U1_1={0.0}
40
41 pressureStep2=LinearSinglePDE(dom) # A [reduced] pressureStep2
42 pressureStep2.setReducedOrderOn() #
43 pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
44 # b free U1={0.0} ?
45
46
47
48 momentumStep3=LinearPDESystem(dom) # A momentumStep3
49 momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
50 # b 1 U1_2={0.0}
51 # b 1 U1_1={0,0}
52 # b 4 U1_1={0.0}
53
54
55 # A deformation
56 # e V700= [grad] V100
57 # e V901= V701+V704 +V101/(X1+small)
58 # e V601=sqrt(2*((V701)^2+(V704)^2+(V101/(X1+small))^2+\
59 # (V702+V703)^2/2))
60 # e U1={abs(V901)/(V601+small)}
61
62
63
64
65 # A stress
66 # e V700= [grad] V100
67 # e V901= V701+V704 +V101/(X1+small)
68 # e V601=sqrt(2*((V701-V901)^2+(V704-V901)^2+(V101/(X1+small)-V901)^2+\
69 # (V702+V703)^2/2))
70 # e V602=((eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small)))
71 # e V704=2*V602*V704-V401
72 # b 3 [integrated] {V704}
73 #
74 # A surface
75 # b 3 [integrated] {1}
76
77
78 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
80 t=dt
81 istep=0
82 En=0
83
84
85 while istep < nstep:
86 istep=istep+1
87 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]
106 # A viscosity
107 gg=grad(U) # V700= [grad] V100
108 vol=gg[0,0]+gg[1,1]+U[0]/(r+small) # V901= V701+V704 +V101/(X1+small)
109 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))
110 # V601=sqrt(2*((V701-V901/3)^2+(V704-V901/3)^2+(V101/(X1+small)-V901/3)^2+(V702+V703)^2/2))
111 m=whereNegative(eta*tau-alpha*p) # eta*V601<alpha*V401
112 eta_d=m*eta+(1.-m)*alpha*p/(tau+small) #
113 # U1={(eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))}
114 # viscosity 1 100 V100 200 V200 300 V300 400 V400
115 # solve
116 print " viscosity =",inf(eta_d),sup(eta_d) # show V801
117 stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
118 # solve momentum equationStep1
119 # momentumStep1 1 100 V100 200 V200 300 V300 400 V400
120 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
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}
123 y=-(nStressn*n_d+nStressTau*t_d))
124 dU2=momentumStep1.getSolution()
125 # solve
126 # dU2-> V100
127 # U -> V200
128 # p -> V500
129
130 # pressureStep2 1 100 V100 200 V200 300 V300 400 V500
131
132 gg2=grad(dU2)
133 div_dU2=gg2[0,0]+gg2[1,1]+dU2[0]/(r+small)
134 pressureStep2.setValue(A=r*teta1*teta2*dt**2*kronecker(dom), # -D_j{teta1*teta2*dt^2}D_jU1
135 D=r*ro, # {ro/Pen}U1
136 Y=-r*(ro*dt*vol+teta1*ro*dt*div_dU2)) # -{ro*dt}D_j{V200}_j-{teta1*ro*dt}D_j{V100}_j
137 # solve
138 dp=pressureStep2.getSolution()
139 # dp -> V100
140 # dU2 -> V200
141 # U -> V300
142 # p -> V600
143 p+=dp # V601=V601+V101
144 # V701=V101 : dp -> V700
145 p=(p+abs(p))/2+small # V601=cut V601
146 # V601=V601+small
147 print " pressure increment :",inf(dp),sup(dp) # show V601
148 print " pressure range :",inf(p),sup(p) # show V601
149 # popp
150 # dU2 -> V100
151 # U -> V200
152 # p -> V500
153 # dp -> V600
154
155 # momentumStep3 1 100 V100 200 V200 300 V300 400 V600
156 A2=Tensor4(0.0,Function(dom))
157 for i in range(dom.getDim()):
158 for j in range(dom.getDim()):
159 A2[i,i,j,j]+=1
160 momentumStep3.setValue(A=Pen*A2, # -D_i{Pen}D_jU1_j
161 D=ro/dt*kronecker(dom), # {ro/dt}U1_i
162 Y=(ro/dt*(dU2+U)-teta2*grad(dp)))
163 # Y=r*(ro/dt*(dU2+U)-teta2*dp*[1.,0]), # {ro/dt}{V200}_i+{ro/dt}{V100}_i
164 # X=r*teta2*dp**kronecker(dom)) # -D_i{teta2*V401}
165 U, U_old=momentumStep3.getSolution(), U
166 # U -> V100
167 # dU2 -> V200
168 # U_old -> V300
169 # p -> V600
170 # dp -> V700
171 # V400=V300
172 # V300=V100
173 # popp
174 # popp
175 # U -> V100
176 # U_old -> V200
177 # p -> V400
178 # dp -> V500
179 # deformation
180 # solve
181 # show ratio
182 # show V101
183 # popp
184
185 print " new U:",inf(U),sup(U) # show v100
186 print " old U",inf(U_old),sup(U_old) # show v200
187 print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200
188 # test0=L1 V100
189 # show test/test0
190 # black
191 # arrow
192
193 vmax=Lsup(U) # vmax1=abs(max V100)
194 # vmax2=abs(min V100)
195 # show vmax1
196 # show vmax2
197
198 # vmax=vmax1
199 # if vmax < vmax2
200 # vmax=vmax2
201 # endif
202 print " max velocity = ",vmax # show vmax
203 len=inf(dom.getSize())
204 dt1=len/vmax # Courant condition
205 dt2=0.5*ro*(len**2)/eta
206 dt=dt1*dt2/(dt1+dt2)
207 En=1/dt
208 t=t+dt
209 print " time , time step ",t,dt # show t
210 # show dt
211
212
213 dom.setX(dom.getX()+U*dt) # mapm 500
214 # V500=V500+V100*dt
215 # V600=V500
216 # mapm 500
217
218 print " volume : ",integrate(r)# vol=Out
219 # show vol
220 # dU -> V100
221 # U -> V200
222 # p -> V400
223 # dp -> V600
224 # clear
225 # prim
226 # saveVTK("u.%d.xml"%(istep),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp)
227 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