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

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-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-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