/[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 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 2994 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 from esys.escript import *
27 from esys.escript.linearPDEs import LinearPDE, TransportPDE
28 from esys.dudley import Rectangle
29 from esys.weipa import saveVTK
30
31 # dom=Rectangle(12,8,l0=1.5)
32 # dom=Rectangle(24,16,l0=1.5)
33 dom=Rectangle(48,32,l0=1.5)
34 # dom=Rectangle(8*48,8*32,l0=1.5)
35 # dom=Rectangle(120,80,l0=1.5)
36 V=Scalar(1.,Function(dom))*[-1.,0]
37 THETA=0.
38 fc=TransportPDE(dom,num_equations=1,theta=THETA)
39 fc.setTolerance(1.e-12)
40 fc.setValue(M=Scalar(1.,Function(dom)),C=V)
41 x=dom.getX()
42 x_0=[0.5,0.5]
43 sigma=0.075
44 u0=1.
45 for i in range(dom.getDim()):
46 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)
47
48 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
49 # f1=0.5
50 # f2=2.
51 # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5
52 # u0=exp(-3*(x[0]-2.)**2)
53 # u0=x[0]
54 u0/=Lsup(u0)
55 c=0
56 saveVTK("u.%s.vtu"%c,u=u0)
57 fc.setInitialSolution(u0)
58
59 t_end=0.6
60 dt=2.49999e-2*0+6.2499999e-02/4
61 dt_out=2.49999e-2*0+6.2499999e-02/4
62 c_stop=1
63 n_out=int(t_end/dt+0.5)
64 print(n_out)
65 t=0.
66 t_out=0
67 c_out=0
68 c=0
69 print(t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom)))
70 while t<t_end and c< c_stop:
71 print("time step t=",t+dt)
72 u=fc.solve(dt)
73 print(t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom)))
74 c+=1
75 t+=dt
76 if t>=t_out+dt_out:
77 c_out,t_out=c_out+1,t_out+dt_out
78 saveVTK("u.%s.vtu"%c_out,u=u)
79 print("write time step ",c,"(t=%s) to file u.%s.vtu"%(t,c_out))
80
81 if True:
82 pde=LinearPDE(dom)
83 pde.setValue(D=1.,C=-THETA*dt*V)
84 pde.setTolerance(1e-12)
85 t=0.
86 t_out=0
87 c_out=0
88 c=0
89 u=u0
90 print(t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom)))
91 while t<t_end and c< c_stop:
92 print("time step t=",t+dt)
93 pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
94 u=pde.getSolution(verbose=True)
95 print(t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom)))
96 c+=1
97 t+=dt
98 if t>=t_out+dt_out:
99 c_out,t_out=c_out+1,t_out+dt_out
100 saveVTK("u2.%s.vtu"%c_out,u=u)
101 print("write time step ",c,"(t=%s) to file u2.%s.vtu"%(t,c_out))

  ViewVC Help
Powered by ViewVC 1.1.26