/[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 1562 - (hide annotations)
Wed May 21 13:04:40 2008 UTC (12 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 3991 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


1 gross 1562 #
2     # AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
3     #
4     #
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 3b p+=1/2(p+dp+abs(p+dp))
9     # step 4 sigma'i+_ij,j=f(v+,p+,...)
10     #
11     #
12     from esys.escript import *
13     from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
14     from esys.finley import Rectangle
15    
16    
17     nel = 20
18     H = 0.5
19     L = 1.0
20    
21     eta = 1.0 # shear viscosity
22     ro = 1.0
23     g = 1.00
24    
25     alpha_w = 1.00
26     alpha = 1.00*1000000.
27     Pen=0.
28     B=100.
29    
30     nstep = 3000
31     dt = 1.
32     small = EPSILON
33     w_step=max(int(nstep/50),1)*0+1
34     toler = 0.001
35     teta1 = 0.5
36     teta2 = 0.5
37     teta3 = 1 # =0 split A; =1 split B
38    
39     # create domain:
40     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    
44     momentumStep1=LinearPDESystem(dom)
45     momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) # fix x0=0 and x1=0
46     face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])
47    
48     pressureStep2=LinearSinglePDE(dom)
49     pressureStep2.setReducedOrderOn()
50     pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
51    
52    
53    
54     momentumStep3=LinearPDESystem(dom)
55     momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])
56     #
57     # initial values:
58     #
59     U=Vector(0.,Solution(dom))
60     p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3
61     # p=ro*g*(H-ReducedSolution(dom).getX()[1])
62     stress=Tensor(0.,Function(dom))
63    
64     t=dt
65     istep=0
66     while istep < nstep:
67     istep=istep+1
68     print "time step :",istep," t = ",t
69     r=Function(dom).getX()[0]
70     r_b=FunctionOnBoundary(dom).getX()[0]
71     print " volume : ",integrate(r)
72     #
73     # step 1:
74     #
75     # calculate normal
76     n_d=dom.getNormal()
77     t_d=matrixmult(numarray.array([[0.,-1.],[1.,0]]),n_d)
78     sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)
79     print " sigma_d =",inf(sigma_d),sup(sigma_d)
80    
81     momentumStep1.setValue(D=ro*kronecker(dom),
82     Y=ro*U+dt*((stress[:,0]-p*teta3*kronecker(dom)[:,0])/r+[0.,-ro*g]),
83     # Y=r*ro*U+dt*r*[0.,-ro*g],
84     X=-dt*(stress-teta3*p*kronecker(dom)),
85     y=sigma_d*face_mask)
86     U_star=momentumStep1.getSolution()
87     #
88     # step 2:
89     #
90     # U2=U+teta1*(U_star-U)
91     U2=U+teta1*U_star
92     gg2=grad(U2)
93     div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r
94    
95     grad_p=grad(p)
96    
97     pressureStep2.setValue(A=r*dt*B*teta1*teta2/ro*dt*kronecker(dom),
98     D=r,
99     Y=-dt*B*r*div_U2,
100     X=-r*B*dt**2/ro*teta1*(1-teta3)*grad_p)
101     dp=pressureStep2.getSolution()
102     #
103     # step 3:
104     #
105     p2=(1-teta3)*p+teta2*dp
106     grad_p2=grad(p2)
107     momentumStep3.setValue(D=r*ro*kronecker(dom),
108     Y=r*(ro*U_star-dt*teta2*grad_p2))
109     U_new=momentumStep3.getSolution()
110     #
111     # update:
112     #
113     p+=dp
114     U=U_new
115     print " U:",inf(U),sup(U)
116     print " P:",inf(p),sup(p)
117    
118    
119     p_pos=clip(p,small)
120     gg=grad(U)
121     vol=gg[0,0]+gg[1,1]+U[0]/r
122     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))
123     m=whereNegative(eta*gamma-alpha*p_pos)
124     eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)
125     print " viscosity =",inf(eta_d),sup(eta_d)
126     stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
127     #
128     # step size control:
129     #
130     len=inf(dom.getSize())
131     dt1=inf(dom.getSize()/(length(U)+small))
132     dt2=inf(0.5*ro*(len**2)/eta_d)
133     dt=dt1*dt2/(dt1+dt2)
134     print " new step size = ",dt
135     #
136     # update geometry
137     #
138     dom.setX(dom.getX()+U*dt)
139     t=t+dt
140     if (istep-1)%w_step==0:saveVTK("u.%d.xml"%((istep-1)/w_step),p=p,eta=eta_d,U=U_star,U_star=U_star,gamma=gamma)
141     if istep == 3: 1/0

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26