← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/axisymm-splitB.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

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