/[escript]/trunk/finley/test/python/tp.py
ViewVC logotype

Contents of /trunk/finley/test/python/tp.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1410 - (show annotations)
Thu Feb 7 04:24:00 2008 UTC (12 years ago) by gross
File MIME type: text/x-python
File size: 2007 byte(s)
a new version of the algebric upwinding a. flux limiter
1 from esys.escript import *
2 from esys.escript.linearPDEs import LinearPDE, TransportPDE
3 from esys.finley import Rectangle
4
5 # 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(8*48,8*32,l0=1.5)
9 # dom=Rectangle(120,80,l0=1.5)
10 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 x=dom.getX()
16 x_0=[0.5,0.5]
17 sigma=0.075
18 u0=1.
19 for i in range(dom.getDim()):
20 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)
21
22 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=f2*clip(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 c=0
30 saveVTK("u.%s.xml"%c,u=u0)
31 fc.setInitialSolution(u0)
32
33 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 t=0.
40 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 print "time step t=",t+dt
46 u=fc.solve(dt)
47 print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
48 c+=1
49 t+=dt
50 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=-THETA*dt*V)
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)*dt*inner(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)

  ViewVC Help
Powered by ViewVC 1.1.26