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

Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 4 months ago) by ksteube
File MIME type: text/x-python
File size: 2772 byte(s)
```Copyright updated in all files

```
 1 ksteube 1811 2 ######################################################## 3 # 4 # Copyright (c) 2003-2008 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2008 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 20 __url__= 21 22 gross 1368 from esys.escript import * 23 gross 1410 from esys.escript.linearPDEs import LinearPDE, TransportPDE 24 gross 1370 from esys.finley import Rectangle 25 26 gross 1410 # dom=Rectangle(12,8,l0=1.5) 27 # dom=Rectangle(24,16,l0=1.5) 28 dom=Rectangle(48,32,l0=1.5) 29 # dom=Rectangle(8*48,8*32,l0=1.5) 30 gross 1407 # dom=Rectangle(120,80,l0=1.5) 31 gross 1410 V=Scalar(1.,Function(dom))*[-1.,0] 32 THETA=0. 33 fc=TransportPDE(dom,num_equations=1,theta=THETA) 34 fc.setTolerance(1.e-12) 35 fc.setValue(M=Scalar(1.,Function(dom)),C=V) 36 gross 1370 x=dom.getX() 37 gross 1410 x_0=[0.5,0.5] 38 sigma=0.075 39 u0=1. 40 gross 1371 for i in range(dom.getDim()): 41 gross 1410 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2) 42 gross 1371 43 gross 1410 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2) 44 # f1=0.5 45 # f2=2. 46 # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5 47 # u0=exp(-3*(x[0]-2.)**2) 48 # u0=x[0] 49 u0/=Lsup(u0) 50 gross 1370 c=0 51 gross 1410 saveVTK("u.%s.xml"%c,u=u0) 52 fc.setInitialSolution(u0) 53 gross 1371 54 gross 1410 t_end=0.6 55 dt=2.49999e-2*0+6.2499999e-02/4 56 dt_out=2.49999e-2*0+6.2499999e-02/4 57 c_stop=1 58 n_out=int(t_end/dt+0.5) 59 print n_out 60 gross 1370 t=0. 61 gross 1410 t_out=0 62 c_out=0 63 c=0 64 print t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom)) 65 while t=t_out+dt_out: 72 c_out,t_out=c_out+1,t_out+dt_out 73 saveVTK("u.%s.xml"%c_out,u=u) 74 print "write time step ",c,"(t=%s) to file u.%s.xml"%(t,c_out) 75 76 if True: 77 pde=LinearPDE(dom) 78 pde.setValue(D=1.,C=-THETA*dt*V) 79 pde.setTolerance(1e-12) 80 t=0. 81 t_out=0 82 c_out=0 83 c=0 84 u=u0 85 print t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom)) 86 while t=t_out+dt_out: 94 c_out,t_out=c_out+1,t_out+dt_out 95 saveVTK("u2.%s.xml"%c_out,u=u) 96 print "write time step ",c,"(t=%s) to file u2.%s.xml"%(t,c_out)