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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 4 months ago) by ksteube
File MIME type: text/x-python
File size: 2772 byte(s)
Copyright updated in all files

1 ksteube 1811
2     ########################################################
3     #
4     # Copyright (c) 2003-2008 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-2008 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__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 gross 1368 from esys.escript import *
23 gross 1410 from esys.escript.linearPDEs import LinearPDE, TransportPDE
24 gross 1370 from esys.finley import Rectangle
25    
26 gross 1410 # dom=Rectangle(12,8,l0=1.5)
27     # dom=Rectangle(24,16,l0=1.5)
28     dom=Rectangle(48,32,l0=1.5)
29     # dom=Rectangle(8*48,8*32,l0=1.5)
30 gross 1407 # dom=Rectangle(120,80,l0=1.5)
31 gross 1410 V=Scalar(1.,Function(dom))*[-1.,0]
32     THETA=0.
33     fc=TransportPDE(dom,num_equations=1,theta=THETA)
34     fc.setTolerance(1.e-12)
35     fc.setValue(M=Scalar(1.,Function(dom)),C=V)
36 gross 1370 x=dom.getX()
37 gross 1410 x_0=[0.5,0.5]
38     sigma=0.075
39     u0=1.
40 gross 1371 for i in range(dom.getDim()):
41 gross 1410 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)
42 gross 1371
43 gross 1410 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
44     # f1=0.5
45     # f2=2.
46     # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5
47     # u0=exp(-3*(x[0]-2.)**2)
48     # u0=x[0]
49     u0/=Lsup(u0)
50 gross 1370 c=0
51 gross 1410 saveVTK("u.%s.xml"%c,u=u0)
52     fc.setInitialSolution(u0)
53 gross 1371
54 gross 1410 t_end=0.6
55     dt=2.49999e-2*0+6.2499999e-02/4
56     dt_out=2.49999e-2*0+6.2499999e-02/4
57     c_stop=1
58     n_out=int(t_end/dt+0.5)
59     print n_out
60 gross 1370 t=0.
61 gross 1410 t_out=0
62     c_out=0
63     c=0
64     print t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom))
65     while t<t_end and c< c_stop:
66 gross 1370 print "time step t=",t+dt
67     u=fc.solve(dt)
68 gross 1410 print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
69 gross 1370 c+=1
70     t+=dt
71 gross 1410 if t>=t_out+dt_out:
72     c_out,t_out=c_out+1,t_out+dt_out
73     saveVTK("u.%s.xml"%c_out,u=u)
74     print "write time step ",c,"(t=%s) to file u.%s.xml"%(t,c_out)
75    
76     if True:
77     pde=LinearPDE(dom)
78     pde.setValue(D=1.,C=-THETA*dt*V)
79     pde.setTolerance(1e-12)
80     t=0.
81     t_out=0
82     c_out=0
83     c=0
84     u=u0
85     print t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom))
86     while t<t_end and c< c_stop:
87     print "time step t=",t+dt
88     pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
89     u=pde.getSolution(verbose=True)
90     print t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom))
91     c+=1
92     t+=dt
93     if t>=t_out+dt_out:
94     c_out,t_out=c_out+1,t_out+dt_out
95     saveVTK("u2.%s.xml"%c_out,u=u)
96     print "write time step ",c,"(t=%s) to file u2.%s.xml"%(t,c_out)

  ViewVC Help
Powered by ViewVC 1.1.26