/[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 1417 - (hide 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 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 gross 1417 w_step=max(int(nstep/50),1)*0+1
10 gross 1392 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 gross 1417 Etau=1000000000.
29 gross 1392
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 gross 1417 momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
35     face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
36    
37 gross 1392 # 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 gross 1417 momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
50 gross 1392 # 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 gross 1417 En=0
83 gross 1392
84    
85     while istep < nstep:
86     istep=istep+1
87     print "time step :",istep," t = ",t
88 gross 1417
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 gross 1392 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 gross 1417 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 gross 1392 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 gross 1417
185 gross 1392 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 gross 1417 En=1/dt
208 gross 1392 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