test/python/axisymm-splitB.py


##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################

from __future__ import print_function, division

__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

#
#   AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION
#
#
#    step 1 rho*(v_star-v) = dt * (sigma'_ij,j-teta3*p,i+f_i)
#    step 2 dp=-dt*B*(v_j,j+teta1*v_star_j,j-dt*teta1*((1-teta3)*p_,jj+teta2*dp_,jj))
#    step 3 rho*(v+-v) = -dt*((1-teta3)*p_,jj+teta2*dp_,jj)
#    step 3b p+=1/2(p+dp+abs(p+dp))
#    step 4 sigma'i+_ij,j=f(v+,p+,...)
#
#
from esys.escript import *
from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem
from esys.finley import Rectangle
from esys.weipa import *


nel      =   20
H        =   0.5
L        =   1.0

eta      =   1.0       # shear viscosity
ro       =   1.0
g        =   1.00

alpha_w   =   1.00
alpha    =   1.00*1000000.
Pen=0.
B=100.

nstep    =   3000
dt       =   1.
small    =   EPSILON
w_step=max(int(nstep/50),1)*0+1
toler    =   0.001
teta1    =    0.5
teta2    =    0.5
teta3    =    1  # =0 split A; =1 split B

# create domain:
dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H)
x=dom.getX()


momentumStep1=LinearPDESystem(dom) 
momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) # fix x0=0 and x1=0
face_mask=whereZero(FunctionOnBoundary(dom).getX()[1])

pressureStep2=LinearSinglePDE(dom) 
pressureStep2.setReducedOrderOn() 
pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))

momentumStep3=LinearPDESystem(dom)
momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])
#
#   initial values:
#
U=Vector(0.,Solution(dom)) 
p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 
p=ro*g*(H-ReducedSolution(dom).getX()[1])
dev_stress=Tensor(0.,Function(dom))

t=dt
istep=0
while istep < nstep:
    istep=istep+1
    print("time step :",istep," t = ",t)
    r=Function(dom).getX()[0]
    r_b=FunctionOnBoundary(dom).getX()[0]
    print("volume : ",integrate(r))
    #
    #  step 1:
    #
    # calculate normal 
    n_d=dom.getNormal()
    t_d=matrixmult(numpy.array([[0.,-1.],[1.,0]]),n_d)
    sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)
    print("sigma_d =",inf(sigma_d),sup(sigma_d))

    momentumStep1.setValue(D=r*ro*kronecker(dom),
                           Y=r*ro*U+dt*r*[0.,-ro*g], 
                           X=-dt*r*(dev_stress-teta3*p*kronecker(dom)), 
                           y=sigma_d*face_mask*r_b)
    U_star=momentumStep1.getSolution()
    saveVTK("u.vtu",u=U_star,u0=U)
    #
    #  step 2:
    #
    # U2=U+teta1*(U_star-U)
    U2=U+teta1*U_star
    gg2=grad(U2)
    div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r

    grad_p=grad(p)

    pressureStep2.setValue(A=r*dt*B*teta1*teta2/ro*dt*kronecker(dom), 
                           D=r,                            
                           Y=-dt*B*r*div_U2,
                           X=-r*B*dt**2/ro*teta1*(1-teta3)*grad_p)
    dp=pressureStep2.getSolution()
    #
    #  step 3:
    #
    p2=(1-teta3)*p+teta2*dp
    grad_p2=grad(p2)
    momentumStep3.setValue(D=r*ro*kronecker(dom),
                           Y=r*(ro*U_star-dt*teta2*grad_p2))
    U_new=momentumStep3.getSolution()
    #
    #   update:
    #
    p+=dp         
    U=U_new
    print("U:",inf(U),sup(U))
    print("P:",inf(p),sup(p)) 


    p_pos=clip(p,small)
    gg=grad(U) 
    vol=gg[0,0]+gg[1,1]+U[0]/r  
    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))
    m=whereNegative(eta*gamma-alpha*p_pos) 
    eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)  
    print("viscosity =",inf(eta_d),sup(eta_d)) 
    dev_stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
    #
    # step size control:
    #
    len=inf(dom.getSize())
    dt1=inf(dom.getSize()/(length(U)+small))
    dt2=inf(0.5*ro*(len**2)/eta_d)
    dt=dt1*dt2/(dt1+dt2)
    print("new step size = ",dt)
    #
    #  update geometry
    #
    dom.setX(dom.getX()+U*dt)
    t=t+dt
    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)
    if istep == 3: 1/0
1	ksteube	1811
2	jfenwick	3981	##############################################################################
3	gross	1562	#
4	jfenwick	6651	# Copyright (c) 2003-2018 by The University of Queensland
5	jfenwick	3981	# http://www.uq.edu.au
6	ksteube	1811	#
7			# Primary Business: Queensland, Australia
8	jfenwick	6112	# Licensed under the Apache License, version 2.0
9			# http://www.apache.org/licenses/LICENSE-2.0
10	ksteube	1811	#
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	6651	__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20	jfenwick	3981	http://www.uq.edu.au
21	ksteube	1811	Primary Business: Queensland, Australia"""
22	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
23			http://www.apache.org/licenses/LICENSE-2.0"""
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=-dtB(v_j,j+teta1v_star_j,j-dtteta1((1-teta3)p_,jj+teta2*dp_,jj))
32			# step 3 rho(v+-v) = -dt((1-teta3)p_,jj+teta2dp_,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(nelL/min(L,H)),int(nelH/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=rog(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3
85	gross	1639	p=rog(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_wt_d-n_d)Penclip(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=rrokronecker(dom),
106			Y=rroU+dtr[0.,-ro*g],
107			X=-dtr(dev_stress-teta3pkronecker(dom)),
108			y=sigma_dface_maskr_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=rdtBteta1teta2/rodtkronecker(dom),
122			D=r,
123			Y=-dtBr*div_U2,
124			X=-rBdt*2/roteta1(1-teta3)grad_p)
125			dp=pressureStep2.getSolution()
126			#
127			# step 3:
128			#
129			p2=(1-teta3)p+teta2dp
130			grad_p2=grad(p2)
131			momentumStep3.setValue(D=rrokronecker(dom),
132			Y=r(roU_star-dtteta2grad_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(etagamma-alphap_pos)
148			eta_d=meta+(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.5ro(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