# Contents of /branches/symbolic_from_3470/dudley/test/python/tp.py

Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 2807 byte(s)
```Fast forward to latest trunk revision 3788.

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2010 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-2010 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.dudley import Rectangle 25 from esys.weipa import saveVTK 26 27 # dom=Rectangle(12,8,l0=1.5) 28 # dom=Rectangle(24,16,l0=1.5) 29 dom=Rectangle(48,32,l0=1.5) 30 # dom=Rectangle(8*48,8*32,l0=1.5) 31 # dom=Rectangle(120,80,l0=1.5) 32 V=Scalar(1.,Function(dom))*[-1.,0] 33 THETA=0. 34 fc=TransportPDE(dom,num_equations=1,theta=THETA) 35 fc.setTolerance(1.e-12) 36 fc.setValue(M=Scalar(1.,Function(dom)),C=V) 37 x=dom.getX() 38 x_0=[0.5,0.5] 39 sigma=0.075 40 u0=1. 41 for i in range(dom.getDim()): 42 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2) 43 44 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2) 45 # f1=0.5 46 # f2=2. 47 # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5 48 # u0=exp(-3*(x[0]-2.)**2) 49 # u0=x[0] 50 u0/=Lsup(u0) 51 c=0 52 saveVTK("u.%s.vtu"%c,u=u0) 53 fc.setInitialSolution(u0) 54 55 t_end=0.6 56 dt=2.49999e-2*0+6.2499999e-02/4 57 dt_out=2.49999e-2*0+6.2499999e-02/4 58 c_stop=1 59 n_out=int(t_end/dt+0.5) 60 print(n_out) 61 t=0. 62 t_out=0 63 c_out=0 64 c=0 65 print(t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom))) 66 while t=t_out+dt_out: 73 c_out,t_out=c_out+1,t_out+dt_out 74 saveVTK("u.%s.vtu"%c_out,u=u) 75 print("write time step ",c,"(t=%s) to file u.%s.vtu"%(t,c_out)) 76 77 if True: 78 pde=LinearPDE(dom) 79 pde.setValue(D=1.,C=-THETA*dt*V) 80 pde.setTolerance(1e-12) 81 t=0. 82 t_out=0 83 c_out=0 84 c=0 85 u=u0 86 print(t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom))) 87 while t=t_out+dt_out: 95 c_out,t_out=c_out+1,t_out+dt_out 96 saveVTK("u2.%s.vtu"%c_out,u=u) 97 print("write time step ",c,"(t=%s) to file u2.%s.vtu"%(t,c_out))