/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 3036 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 ksteube 1811
2 jfenwick 3981 ##############################################################################
3 ksteube 1811 #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1811 #
7     # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1811 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 ksteube 1811
17 sshaw 5706 from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 ksteube 1811 Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1811
26 gross 1368 from esys.escript import *
27 gross 1410 from esys.escript.linearPDEs import LinearPDE, TransportPDE
28 gross 1370 from esys.finley import Rectangle
29 caltinay 3346 from esys.weipa import saveVTK
30 gross 1370
31 gross 1410 # dom=Rectangle(12,8,l0=1.5)
32     # dom=Rectangle(24,16,l0=1.5)
33     dom=Rectangle(48,32,l0=1.5)
34 gross 4185 saveDataCSV("t.csv",x=dom.getX(), rho=length(dom.getX()))
35     1/0
36 gross 1410 # dom=Rectangle(8*48,8*32,l0=1.5)
37 gross 1407 # dom=Rectangle(120,80,l0=1.5)
38 gross 1410 V=Scalar(1.,Function(dom))*[-1.,0]
39     THETA=0.
40     fc=TransportPDE(dom,num_equations=1,theta=THETA)
41     fc.setTolerance(1.e-12)
42     fc.setValue(M=Scalar(1.,Function(dom)),C=V)
43 gross 1370 x=dom.getX()
44 gross 1410 x_0=[0.5,0.5]
45     sigma=0.075
46     u0=1.
47 gross 1371 for i in range(dom.getDim()):
48 sshaw 4576 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)
49 gross 1371
50 gross 1410 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
51     # f1=0.5
52     # f2=2.
53     # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5
54     # u0=exp(-3*(x[0]-2.)**2)
55     # u0=x[0]
56     u0/=Lsup(u0)
57 gross 1370 c=0
58 caltinay 2534 saveVTK("u.%s.vtu"%c,u=u0)
59 gross 1410 fc.setInitialSolution(u0)
60 gross 1371
61 gross 1410 t_end=0.6
62     dt=2.49999e-2*0+6.2499999e-02/4
63     dt_out=2.49999e-2*0+6.2499999e-02/4
64     c_stop=1
65     n_out=int(t_end/dt+0.5)
66 jfenwick 3772 print(n_out)
67 gross 1370 t=0.
68 gross 1410 t_out=0
69     c_out=0
70     c=0
71 jfenwick 3772 print(t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom)))
72 gross 1410 while t<t_end and c< c_stop:
73 sshaw 4576 print("time step t=",t+dt)
74     u=fc.solve(dt)
75 jfenwick 3772 print(t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom)))
76 gross 1370 c+=1
77     t+=dt
78 gross 1410 if t>=t_out+dt_out:
79     c_out,t_out=c_out+1,t_out+dt_out
80 caltinay 2534 saveVTK("u.%s.vtu"%c_out,u=u)
81 jfenwick 3772 print("write time step ",c,"(t=%s) to file u.%s.vtu"%(t,c_out))
82 gross 1410
83     if True:
84     pde=LinearPDE(dom)
85     pde.setValue(D=1.,C=-THETA*dt*V)
86     pde.setTolerance(1e-12)
87     t=0.
88     t_out=0
89     c_out=0
90     c=0
91     u=u0
92 jfenwick 3772 print(t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom)))
93 gross 1410 while t<t_end and c< c_stop:
94 sshaw 4576 print("time step t=",t+dt)
95 gross 1410 pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
96     u=pde.getSolution(verbose=True)
97 jfenwick 3772 print(t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom)))
98 gross 1410 c+=1
99     t+=dt
100     if t>=t_out+dt_out:
101     c_out,t_out=c_out+1,t_out+dt_out
102 caltinay 2534 saveVTK("u2.%s.vtu"%c_out,u=u)
103 jfenwick 3772 print("write time step ",c,"(t=%s) to file u2.%s.vtu"%(t,c_out))

  ViewVC Help
Powered by ViewVC 1.1.26