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

Revision 2534 - (show annotations)
Thu Jul 16 06:49:19 2009 UTC (10 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 2767 byte(s)
Changed examples, tests and tutorials to save VTK files as .vtu instead .xml.
Visit doesn't know what to do with xml's and vtu is the proper extension
anyway.

 1 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 from esys.escript import * 23 from esys.escript.linearPDEs import LinearPDE, TransportPDE 24 from esys.finley import Rectangle 25 26 # 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 # dom=Rectangle(120,80,l0=1.5) 31 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 x=dom.getX() 37 x_0=[0.5,0.5] 38 sigma=0.075 39 u0=1. 40 for i in range(dom.getDim()): 41 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2) 42 43 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 c=0 51 saveVTK("u.%s.vtu"%c,u=u0) 52 fc.setInitialSolution(u0) 53 54 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 t=0. 61 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.vtu"%c_out,u=u) 74 print "write time step ",c,"(t=%s) to file u.%s.vtu"%(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.vtu"%c_out,u=u) 96 print "write time step ",c,"(t=%s) to file u2.%s.vtu"%(t,c_out)