/[escript]/trunk/finley/test/python/tp.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1370 by gross, Wed Jan 2 09:21:43 2008 UTC revision 1392 by gross, Mon Jan 21 06:20:54 2008 UTC
# Line 56  class TransportPDE(object): Line 56  class TransportPDE(object):
56       def setInitialSolution(self,u):       def setInitialSolution(self,u):
57           self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))           self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
58       def solve(self,dt):       def solve(self,dt):
59         return self.__transport_problem.solve(self.__source,dt,{})         return self.__transport_problem.solve(self.__source,dt,{"verbose" : True })
60  from esys.finley import Rectangle  from esys.finley import Rectangle
61    
62  dom=Rectangle(10,10)  dom=Rectangle(40,20,l0=2)
63  fc=TransportPDE(dom,num_equations=1,theta=0.0,dt_max=0.)  fc=TransportPDE(dom,num_equations=1,theta=1.0,dt_max=2.5e-2/10)
64  fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0.])  fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])
65  x=dom.getX()  x=dom.getX()
66  u=whereNonPositive(abs(x[0]-0.3)-0.1)*whereNonPositive(abs(x[1]-0.3)-0.1)  x_0=[0.3,0.3]
67    sigma=0.08
68    u=1.
69    for i in range(dom.getDim()):
70        u=u*exp(-(x[i]-x_0[i])**2/sigma**2)
71    u/=Lsup(u)
72    
73    u=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
74  c=0  c=0
75  saveVTK("u.%s.xml"%c,u=u)  saveVTK("u.%s.xml"%c,u=u)
76  fc.setInitialSolution(u)  fc.setInitialSolution(u)
77  dt=0.301  
78    dt=2.5e-2
79  t=0.  t=0.
80  while t<dt:  while t<25*dt:
81      print "time step t=",t+dt        print "time step t=",t+dt  
82      u=fc.solve(dt)        u=fc.solve(dt)  
83      print "range u",inf(u),sup(u)      print "range u",inf(u),sup(u),integrate(u,Function(dom))
84      c+=1      c+=1
85      saveVTK("u.%s.xml"%c,u=u)      saveVTK("u.%s.xml"%c,u=u)
86      t+=dt      t+=dt

Legend:
Removed from v.1370  
changed lines
  Added in v.1392

  ViewVC Help
Powered by ViewVC 1.1.26