/[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 1552 - (hide 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 gross 1392 #
2     # AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
3     #
4 gross 1552 #
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 gross 1392 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 gross 1417 w_step=max(int(nstep/50),1)*0+1
18 gross 1392 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 gross 1552 teta3 = 0 # =0 split A; =1 split B
37 gross 1417 Etau=1000000000.
38 gross 1392
39 gross 1552 # create domain:
40 gross 1392 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 gross 1552
44 gross 1392 momentumStep1=LinearPDESystem(dom) # A momentumStep1
45 gross 1417 momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
46     face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
47    
48 gross 1392 # 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 gross 1417 momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.])
61 gross 1392 # 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 gross 1417 En=0
94 gross 1392
95    
96     while istep < nstep:
97     istep=istep+1
98     print "time step :",istep," t = ",t
99 gross 1417
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 gross 1392 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 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}
134     y=-(nStressn*n_d+nStressTau*t_d))
135 gross 1392 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 gross 1417
196 gross 1392 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 gross 1417 En=1/dt
219 gross 1392 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