/[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 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4795 byte(s)
Round 1 of copyright fixes
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 #
24 # AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
25 #
26 #
27 # step 1 rho*(v_star-v) = dt * (sigma'_ij,j-teta3*p,i+f_i)
28 # step 2 dp=-dt*B*(v_j,j+teta1*v_star_j,j-dt*teta1*((1-teta3)*p_,jj+teta2*dp_,jj))
29 # step 3 rho*(v+-v) = -dt*((1-teta3)*p_,jj+teta2*dp_,jj)
30 # step 3b p+=1/2(p+dp+abs(p+dp))
31 # step 4 sigma'i+_ij,j=f(v+,p+,...)
32 #
33 #
34 from esys.escript import *
35 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
36 from esys.finley import Rectangle
37 from esys.weipa import *
38
39
40 nel = 20
41 H = 0.5
42 L = 1.0
43
44 eta = 1.0 # shear viscosity
45 ro = 1.0
46 g = 1.00
47
48 alpha_w = 1.00
49 alpha = 1.00*1000000.
50 Pen=0.
51 B=100.
52
53 nstep = 3000
54 dt = 1.
55 small = EPSILON
56 w_step=max(int(nstep/50),1)*0+1
57 toler = 0.001
58 teta1 = 0.5
59 teta2 = 0.5
60 teta3 = 1 # =0 split A; =1 split B
61
62 # create domain:
63 dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)
64 x=dom.getX()
65
66
67 momentumStep1=LinearPDESystem(dom)
68 momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) # fix x0=0 and x1=0
69 face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
70
71 pressureStep2=LinearSinglePDE(dom)
72 pressureStep2.setReducedOrderOn()
73 pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
74
75 momentumStep3=LinearPDESystem(dom)
76 momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])
77 #
78 # initial values:
79 #
80 U=Vector(0.,Solution(dom))
81 p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3
82 p=ro*g*(H-ReducedSolution(dom).getX()[1])
83 dev_stress=Tensor(0.,Function(dom))
84
85 t=dt
86 istep=0
87 while istep < nstep:
88 istep=istep+1
89 print("time step :",istep," t = ",t)
90 r=Function(dom).getX()[0]
91 r_b=FunctionOnBoundary(dom).getX()[0]
92 print(" volume : ",integrate(r))
93 #
94 # step 1:
95 #
96 # calculate normal
97 n_d=dom.getNormal()
98 t_d=matrixmult(numpy.array([[0.,-1.],[1.,0]]),n_d)
99 sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)
100 print(" sigma_d =",inf(sigma_d),sup(sigma_d))
101
102 momentumStep1.setValue(D=r*ro*kronecker(dom),
103 Y=r*ro*U+dt*r*[0.,-ro*g],
104 X=-dt*r*(dev_stress-teta3*p*kronecker(dom)),
105 y=sigma_d*face_mask*r_b)
106 U_star=momentumStep1.getSolution()
107 saveVTK("u.vtu",u=U_star,u0=U)
108 #
109 # step 2:
110 #
111 # U2=U+teta1*(U_star-U)
112 U2=U+teta1*U_star
113 gg2=grad(U2)
114 div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r
115
116 grad_p=grad(p)
117
118 pressureStep2.setValue(A=r*dt*B*teta1*teta2/ro*dt*kronecker(dom),
119 D=r,
120 Y=-dt*B*r*div_U2,
121 X=-r*B*dt**2/ro*teta1*(1-teta3)*grad_p)
122 dp=pressureStep2.getSolution()
123 #
124 # step 3:
125 #
126 p2=(1-teta3)*p+teta2*dp
127 grad_p2=grad(p2)
128 momentumStep3.setValue(D=r*ro*kronecker(dom),
129 Y=r*(ro*U_star-dt*teta2*grad_p2))
130 U_new=momentumStep3.getSolution()
131 #
132 # update:
133 #
134 p+=dp
135 U=U_new
136 print(" U:",inf(U),sup(U))
137 print(" P:",inf(p),sup(p))
138
139
140 p_pos=clip(p,small)
141 gg=grad(U)
142 vol=gg[0,0]+gg[1,1]+U[0]/r
143 gamma=sqrt(2*((gg[0,0]-vol/3)**2+(gg[1,1]-vol/3)**2+(U[0]/r-vol/3)**2+(gg[1,0]+gg[0,1])**2/2))
144 m=whereNegative(eta*gamma-alpha*p_pos)
145 eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)
146 print(" viscosity =",inf(eta_d),sup(eta_d))
147 dev_stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
148 #
149 # step size control:
150 #
151 len=inf(dom.getSize())
152 dt1=inf(dom.getSize()/(length(U)+small))
153 dt2=inf(0.5*ro*(len**2)/eta_d)
154 dt=dt1*dt2/(dt1+dt2)
155 print(" new step size = ",dt)
156 #
157 # update geometry
158 #
159 dom.setX(dom.getX()+U*dt)
160 t=t+dt
161 if (istep-1)%w_step==0:saveVTK("u.%d.vtu"%((istep-1)/w_step),p=p,eta=eta_d,U=U_star,U_star=U_star,gamma=gamma)
162 if istep == 3: 1/0

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26