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

Revision 1410 - (hide annotations)
Thu Feb 7 04:24:00 2008 UTC (13 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 2007 byte(s)
```a new version of the algebric upwinding a. flux limiter
```
 1 gross 1368 from esys.escript import * 2 gross 1410 from esys.escript.linearPDEs import LinearPDE, TransportPDE 3 gross 1370 from esys.finley import Rectangle 4 5 gross 1410 # dom=Rectangle(12,8,l0=1.5) 6 # dom=Rectangle(24,16,l0=1.5) 7 dom=Rectangle(48,32,l0=1.5) 8 # dom=Rectangle(8*48,8*32,l0=1.5) 9 gross 1407 # dom=Rectangle(120,80,l0=1.5) 10 gross 1410 V=Scalar(1.,Function(dom))*[-1.,0] 11 THETA=0. 12 fc=TransportPDE(dom,num_equations=1,theta=THETA) 13 fc.setTolerance(1.e-12) 14 fc.setValue(M=Scalar(1.,Function(dom)),C=V) 15 gross 1370 x=dom.getX() 16 gross 1410 x_0=[0.5,0.5] 17 sigma=0.075 18 u0=1. 19 gross 1371 for i in range(dom.getDim()): 20 gross 1410 u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2) 21 gross 1371 22 gross 1410 u0=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2) 23 # f1=0.5 24 # f2=2. 25 # u0=f2*clip(x[0]-0.5,0.)-clip(0.5-x[0],0.)*f1+f1*0.5 26 # u0=exp(-3*(x[0]-2.)**2) 27 # u0=x[0] 28 u0/=Lsup(u0) 29 gross 1370 c=0 30 gross 1410 saveVTK("u.%s.xml"%c,u=u0) 31 fc.setInitialSolution(u0) 32 gross 1371 33 gross 1410 t_end=0.6 34 dt=2.49999e-2*0+6.2499999e-02/4 35 dt_out=2.49999e-2*0+6.2499999e-02/4 36 c_stop=1 37 n_out=int(t_end/dt+0.5) 38 print n_out 39 gross 1370 t=0. 40 gross 1410 t_out=0 41 c_out=0 42 c=0 43 print t,": range u",inf(u0),sup(u0),integrate(u0,Function(dom)) 44 while t=t_out+dt_out: 51 c_out,t_out=c_out+1,t_out+dt_out 52 saveVTK("u.%s.xml"%c_out,u=u) 53 print "write time step ",c,"(t=%s) to file u.%s.xml"%(t,c_out) 54 55 if True: 56 pde=LinearPDE(dom) 57 pde.setValue(D=1.,C=-THETA*dt*V) 58 pde.setTolerance(1e-12) 59 t=0. 60 t_out=0 61 c_out=0 62 c=0 63 u=u0 64 print t,": range u2",inf(u0),sup(u0),integrate(u0,Function(dom)) 65 while t=t_out+dt_out: 73 c_out,t_out=c_out+1,t_out+dt_out 74 saveVTK("u2.%s.xml"%c_out,u=u) 75 print "write time step ",c,"(t=%s) to file u2.%s.xml"%(t,c_out)