/[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 1552 - (show annotations)
Thu May 8 08:52:41 2008 UTC (12 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 8604 byte(s)
some changes to make the implementatiopn of a upwind MPI version easier
1 #
2 # 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)
6 # 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)
8 # step 4 p+=1/2(p+dp+abs(p+dp))
9 # step 4 sigma'i+_ij,j=f(v+,p+,...)
10 #
11 #
12 from esys.escript import *
13 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
14 from esys.finley import Rectangle
15 iter = 0
16 nstep = 3000
17 w_step=max(int(nstep/50),1)*0+1
18 mstep = 5
19 H = 0.5
20 eta = 1.0 # shear viscosity
21 ro = 1.0
22 g = 1.00
23 tauy0 = 0.9*ro*g*H/sqrt(3) #0.11/(3*sqrt(3))
24 Pen = 100
25 vtop = 0.01
26 dt = 10**(-10)
27 small = EPSILON
28 alpha = 1.00
29 nel = 20
30 # len = sqrt(0.5/nel)
31 toler = 0.001
32 alphaw = 1.00
33 L = 1.0
34 teta1 = 0.5
35 teta2 = 0.5
36 teta3 = 0 # =0 split A; =1 split B
37 Etau=1000000000.
38
39 # create domain:
40 dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)
41 x=dom.getX()
42
43
44 momentumStep1=LinearPDESystem(dom) # A momentumStep1
45 momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
46 face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
47
48 # b 1 U1_2={0.0}
49 # b 1 U1_1={0,0}
50 # b 4 U1_1={0.0}
51
52 pressureStep2=LinearSinglePDE(dom) # A [reduced] pressureStep2
53 pressureStep2.setReducedOrderOn() #
54 pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
55 # b free U1={0.0} ?
56
57
58
59 momentumStep3=LinearPDESystem(dom) # A momentumStep3
60 momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
61 # b 1 U1_2={0.0}
62 # b 1 U1_1={0,0}
63 # b 4 U1_1={0.0}
64
65
66 # A deformation
67 # e V700= [grad] V100
68 # e V901= V701+V704 +V101/(X1+small)
69 # e V601=sqrt(2*((V701)^2+(V704)^2+(V101/(X1+small))^2+\
70 # (V702+V703)^2/2))
71 # e U1={abs(V901)/(V601+small)}
72
73
74
75
76 # A stress
77 # e V700= [grad] V100
78 # e V901= V701+V704 +V101/(X1+small)
79 # e V601=sqrt(2*((V701-V901)^2+(V704-V901)^2+(V101/(X1+small)-V901)^2+\
80 # (V702+V703)^2/2))
81 # e V602=((eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small)))
82 # e V704=2*V602*V704-V401
83 # b 3 [integrated] {V704}
84 #
85 # A surface
86 # b 3 [integrated] {1}
87
88
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
91 t=dt
92 istep=0
93 En=0
94
95
96 while istep < nstep:
97 istep=istep+1
98 print "time step :",istep," t = ",t
99
100 n_d=dom.getNormal()
101 t_d=matrixmult(numarray.array([[0.,1.],[-1.,0]]),n_d)
102
103 s=-En*inner(U,n_d)*face_mask
104 nStressn=s*wherePositive(s)
105 u_boundary=Etau*inner(U,t_d)*face_mask
106 m=whereNegative(u_boundary-nStressn)
107 nStressTau=nStressn*(1.-m)+u_boundary*m
108 print nStressn
109 # print nStressn*n_d+nStressTau*t_d
110 if istep == 20:
111 # print nStressTau
112 saveVTK("stress.vtu",s=(nStressn*n_d+nStressTau*t_d))
113 1/0
114 # print nStressTau
115
116 r=dom.getX()[0]
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)
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))
121 # V601=sqrt(2*((V701-V901/3)^2+(V704-V901/3)^2+(V101/(X1+small)-V901/3)^2+(V702+V703)^2/2))
122 m=whereNegative(eta*tau-alpha*p) # eta*V601<alpha*V401
123 eta_d=m*eta+(1.-m)*alpha*p/(tau+small) #
124 # U1={(eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))}
125 # viscosity 1 100 V100 200 V200 300 V300 400 V400
126 # solve
127 print " viscosity =",inf(eta_d),sup(eta_d) # show V801
128 stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
129 # solve momentum equationStep1
130 # momentumStep1 1 100 V100 200 V200 300 V300 400 V400
131 momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i
132 Y=stress[:,0]/(r+small)+[0.,-ro*g], # -{0.0,ro*g}_i
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}
134 y=-(nStressn*n_d+nStressTau*t_d))
135 dU2=momentumStep1.getSolution()
136 # solve
137 # dU2-> V100
138 # U -> V200
139 # p -> V500
140
141 # pressureStep2 1 100 V100 200 V200 300 V300 400 V500
142
143 gg2=grad(dU2)
144 div_dU2=gg2[0,0]+gg2[1,1]+dU2[0]/(r+small)
145 pressureStep2.setValue(A=r*teta1*teta2*dt**2*kronecker(dom), # -D_j{teta1*teta2*dt^2}D_jU1
146 D=r*ro, # {ro/Pen}U1
147 Y=-r*(ro*dt*vol+teta1*ro*dt*div_dU2)) # -{ro*dt}D_j{V200}_j-{teta1*ro*dt}D_j{V100}_j
148 # solve
149 dp=pressureStep2.getSolution()
150 # dp -> V100
151 # dU2 -> V200
152 # U -> V300
153 # p -> V600
154 p+=dp # V601=V601+V101
155 # V701=V101 : dp -> V700
156 p=(p+abs(p))/2+small # V601=cut V601
157 # V601=V601+small
158 print " pressure increment :",inf(dp),sup(dp) # show V601
159 print " pressure range :",inf(p),sup(p) # show V601
160 # popp
161 # dU2 -> V100
162 # U -> V200
163 # p -> V500
164 # dp -> V600
165
166 # momentumStep3 1 100 V100 200 V200 300 V300 400 V600
167 A2=Tensor4(0.0,Function(dom))
168 for i in range(dom.getDim()):
169 for j in range(dom.getDim()):
170 A2[i,i,j,j]+=1
171 momentumStep3.setValue(A=Pen*A2, # -D_i{Pen}D_jU1_j
172 D=ro/dt*kronecker(dom), # {ro/dt}U1_i
173 Y=(ro/dt*(dU2+U)-teta2*grad(dp)))
174 # Y=r*(ro/dt*(dU2+U)-teta2*dp*[1.,0]), # {ro/dt}{V200}_i+{ro/dt}{V100}_i
175 # X=r*teta2*dp**kronecker(dom)) # -D_i{teta2*V401}
176 U, U_old=momentumStep3.getSolution(), U
177 # U -> V100
178 # dU2 -> V200
179 # U_old -> V300
180 # p -> V600
181 # dp -> V700
182 # V400=V300
183 # V300=V100
184 # popp
185 # popp
186 # U -> V100
187 # U_old -> V200
188 # p -> V400
189 # dp -> V500
190 # deformation
191 # solve
192 # show ratio
193 # show V101
194 # popp
195
196 print " new U:",inf(U),sup(U) # show v100
197 print " old U",inf(U_old),sup(U_old) # show v200
198 print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200
199 # test0=L1 V100
200 # show test/test0
201 # black
202 # arrow
203
204 vmax=Lsup(U) # vmax1=abs(max V100)
205 # vmax2=abs(min V100)
206 # show vmax1
207 # show vmax2
208
209 # vmax=vmax1
210 # if vmax < vmax2
211 # vmax=vmax2
212 # endif
213 print " max velocity = ",vmax # show vmax
214 len=inf(dom.getSize())
215 dt1=len/vmax # Courant condition
216 dt2=0.5*ro*(len**2)/eta
217 dt=dt1*dt2/(dt1+dt2)
218 En=1/dt
219 t=t+dt
220 print " time , time step ",t,dt # show t
221 # show dt
222
223
224 dom.setX(dom.getX()+U*dt) # mapm 500
225 # V500=V500+V100*dt
226 # V600=V500
227 # mapm 500
228
229 print " volume : ",integrate(r)# vol=Out
230 # show vol
231 # dU -> V100
232 # U -> V200
233 # p -> V400
234 # dp -> V600
235 # clear
236 # prim
237 # saveVTK("u.%d.xml"%(istep),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp)
238 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