/[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 1388 by trankine, Fri Jan 11 07:45:58 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(20,20)  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  x_0=[0.3,0.3]  x_0=[0.3,0.3]
# Line 70  for i in range(dom.getDim()): Line 70  for i in range(dom.getDim()):
70      u=u*exp(-(x[i]-x_0[i])**2/sigma**2)      u=u*exp(-(x[i]-x_0[i])**2/sigma**2)
71  u/=Lsup(u)  u/=Lsup(u)
72    
73  # u=whereNonPositive(abs(x[0]-0.3)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)  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    
78  dt=2.5e-2  dt=2.5e-2
79  t=0.  t=0.
80  while t<15*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),integrate(u,Function(dom))      print "range u",inf(u),sup(u),integrate(u,Function(dom))
Legend:
Removed from v.1388  
changed lines
  Added in v.1392
  ViewVC Help
Powered by ViewVC 1.1.26