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

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 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 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 24 __url__= 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_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_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))