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

Revision 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 1 month ago) by sshaw
File MIME type: text/x-python
File size: 2931 byte(s)
```python3ified things, replaced mixed whitespace and xrange calls
```
 1 2 ############################################################################## 3 # 4 # Copyright (c) 2003-2013 by University of Queensland 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 # Development since 2012 by School of Earth Sciences 13 # 14 ############################################################################## 15 16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland 17 http://www.uq.edu.au 18 Primary Business: Queensland, Australia""" 19 __license__="""Licensed under the Open Software License version 3.0 20 21 __url__= 22 23 from esys.escript import * 24 from esys.escript.linearPDEs import LinearPDE, TransportPDE 25 from esys.finley import Rectangle 26 from esys.weipa import saveVTK 27 28 # dom=Rectangle(12,8,l0=1.5) 29 # dom=Rectangle(24,16,l0=1.5) 30 dom=Rectangle(48,32,l0=1.5) 31 saveDataCSV("t.csv",x=dom.getX(), rho=length(dom.getX())) 32 1/0 33 # dom=Rectangle(8*48,8*32,l0=1.5) 34 # dom=Rectangle(120,80,l0=1.5) 35 V=Scalar(1.,Function(dom))*[-1.,0] 36 THETA=0. 37 fc=TransportPDE(dom,num_equations=1,theta=THETA) 38 fc.setTolerance(1.e-12) 39 fc.setValue(M=Scalar(1.,Function(dom)),C=V) 40 x=dom.getX() 41 x_0=[0.5,0.5] 42 sigma=0.075 43 u0=1. 44 for i in range(dom.getDim()): 45 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2) 46 47 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2) 48 # f1=0.5 49 # f2=2. 50 # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5 51 # u0=exp(-3*(x[0]-2.)**2) 52 # u0=x[0] 53 u0/=Lsup(u0) 54 c=0 55 saveVTK("u.%s.vtu"%c,u=u0) 56 fc.setInitialSolution(u0) 57 58 t_end=0.6 59 dt=2.49999e-2*0+6.2499999e-02/4 60 dt_out=2.49999e-2*0+6.2499999e-02/4 61 c_stop=1 62 n_out=int(t_end/dt+0.5) 63 print(n_out) 64 t=0. 65 t_out=0 66 c_out=0 67 c=0 68 print(t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom))) 69 while t=t_out+dt_out: 76 c_out,t_out=c_out+1,t_out+dt_out 77 saveVTK("u.%s.vtu"%c_out,u=u) 78 print("write time step ",c,"(t=%s) to file u.%s.vtu"%(t,c_out)) 79 80 if True: 81 pde=LinearPDE(dom) 82 pde.setValue(D=1.,C=-THETA*dt*V) 83 pde.setTolerance(1e-12) 84 t=0. 85 t_out=0 86 c_out=0 87 c=0 88 u=u0 89 print(t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom))) 90 while t=t_out+dt_out: 98 c_out,t_out=c_out+1,t_out+dt_out 99 saveVTK("u2.%s.vtu"%c_out,u=u) 100 print("write time step ",c,"(t=%s) to file u2.%s.vtu"%(t,c_out))