test/python/tp.py

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE, TransportPDE
from esys.finley import Rectangle

# dom=Rectangle(12,8,l0=1.5)
# dom=Rectangle(24,16,l0=1.5)
dom=Rectangle(48,32,l0=1.5)
# dom=Rectangle(8*48,8*32,l0=1.5)
# dom=Rectangle(120,80,l0=1.5)
V=Scalar(1.,Function(dom))*[-1.,0]
THETA=0.
fc=TransportPDE(dom,num_equations=1,theta=THETA)
fc.setTolerance(1.e-12)
fc.setValue(M=Scalar(1.,Function(dom)),C=V)
x=dom.getX()
x_0=[0.5,0.5]
sigma=0.075
u0=1.
for i in range(dom.getDim()):
    u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)

u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
# f1=0.5
# f2=2.
# u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5
# u0=exp(-3*(x[0]-2.)**2)
# u0=x[0]
u0/=Lsup(u0)
c=0
saveVTK("u.%s.xml"%c,u=u0)
fc.setInitialSolution(u0)

t_end=0.6
dt=2.49999e-2*0+6.2499999e-02/4
dt_out=2.49999e-2*0+6.2499999e-02/4
c_stop=1
n_out=int(t_end/dt+0.5)
print n_out
t=0.
t_out=0
c_out=0
c=0
print t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom))
while t<t_end and c< c_stop:
    print "time step t=",t+dt   
    u=fc.solve(dt)  
    print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
    c+=1
    t+=dt
    if t>=t_out+dt_out:
         c_out,t_out=c_out+1,t_out+dt_out
         saveVTK("u.%s.xml"%c_out,u=u)
         print "write time step ",c,"(t=%s) to file u.%s.xml"%(t,c_out)

if True:
   pde=LinearPDE(dom)
   pde.setValue(D=1.,C=-THETA*dt*V)
   pde.setTolerance(1e-12)
   t=0.
   t_out=0
   c_out=0
   c=0
   u=u0
   print t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom))
   while t<t_end and c< c_stop:
       print "time step t=",t+dt    
       pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
       u=pde.getSolution(verbose=True)
       print t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom))
       c+=1
       t+=dt
       if t>=t_out+dt_out:
         c_out,t_out=c_out+1,t_out+dt_out
         saveVTK("u2.%s.xml"%c_out,u=u)
         print "write time step ",c,"(t=%s) to file u2.%s.xml"%(t,c_out)
1	gross	1368	from esys.escript import *
2	gross	1410	from esys.escript.linearPDEs import LinearPDE, TransportPDE
3	gross	1370	from esys.finley import Rectangle
4
5	gross	1410	# dom=Rectangle(12,8,l0=1.5)
6			# dom=Rectangle(24,16,l0=1.5)
7			dom=Rectangle(48,32,l0=1.5)
8			# dom=Rectangle(848,832,l0=1.5)
9	gross	1407	# dom=Rectangle(120,80,l0=1.5)
10	gross	1410	V=Scalar(1.,Function(dom))*[-1.,0]
11			THETA=0.
12			fc=TransportPDE(dom,num_equations=1,theta=THETA)
13			fc.setTolerance(1.e-12)
14			fc.setValue(M=Scalar(1.,Function(dom)),C=V)
15	gross	1370	x=dom.getX()
16	gross	1410	x_0=[0.5,0.5]
17			sigma=0.075
18			u0=1.
19	gross	1371	for i in range(dom.getDim()):
20	gross	1410	u0=u0exp(-(x[i]-x_0[i])2/sigma*2)
21	gross	1371
22	gross	1410	u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
23			# f1=0.5
24			# f2=2.
25			# u0=f2clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)f1+f1*0.5
26			# u0=exp(-3(x[0]-2.)*2)
27			# u0=x[0]
28			u0/=Lsup(u0)
29	gross	1370	c=0
30	gross	1410	saveVTK("u.%s.xml"%c,u=u0)
31			fc.setInitialSolution(u0)
32	gross	1371
33	gross	1410	t_end=0.6
34			dt=2.49999e-2*0+6.2499999e-02/4
35			dt_out=2.49999e-2*0+6.2499999e-02/4
36			c_stop=1
37			n_out=int(t_end/dt+0.5)
38			print n_out
39	gross	1370	t=0.
40	gross	1410	t_out=0
41			c_out=0
42			c=0
43			print t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom))
44			while t<t_end and c< c_stop:
45	gross	1370	print "time step t=",t+dt
46			u=fc.solve(dt)
47	gross	1410	print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
48	gross	1370	c+=1
49			t+=dt
50	gross	1410	if t>=t_out+dt_out:
51			c_out,t_out=c_out+1,t_out+dt_out
52			saveVTK("u.%s.xml"%c_out,u=u)
53			print "write time step ",c,"(t=%s) to file u.%s.xml"%(t,c_out)
54
55			if True:
56			pde=LinearPDE(dom)
57			pde.setValue(D=1.,C=-THETAdtV)
58			pde.setTolerance(1e-12)
59			t=0.
60			t_out=0
61			c_out=0
62			c=0
63			u=u0
64			print t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom))
65			while t<t_end and c< c_stop:
66			print "time step t=",t+dt
67			pde.setValue(Y=u+(1.-THETA)dtinner(V,grad(u)))
68			u=pde.getSolution(verbose=True)
69			print t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom))
70			c+=1
71			t+=dt
72			if t>=t_out+dt_out:
73			c_out,t_out=c_out+1,t_out+dt_out
74			saveVTK("u2.%s.xml"%c_out,u=u)
75			print "write time step ",c,"(t=%s) to file u2.%s.xml"%(t,c_out)