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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 2807 byte(s)
dudley, pasowrap, pycad

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

  ViewVC Help
Powered by ViewVC 1.1.26