/[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 5863 - (hide annotations)
Wed Jan 13 02:25:48 2016 UTC (3 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 4914 byte(s)
Copyright dates updated.
\version for doxygen to read

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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26