/[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 1392 by gross, Mon Jan 21 06:20:54 2008 UTC revision 2534 by caltinay, Thu Jul 16 06:49:19 2009 UTC
# Line 1  Line 1 
 from esys.escript import *  
1    
2  class TransportPDE(object):  ########################################################
3       def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):  #
4          self.__domain=domain  # Copyright (c) 2003-2008 by University of Queensland
5          self.__num_equations=num_equations  # Earth Systems Science Computational Center (ESSCC)
6          self.__theta=theta  # http://www.uq.edu.au/esscc
7          self.__dt_max=dt_max  #
8          self.__transport_problem=None  # Primary Business: Queensland, Australia
9      self.__trace=trace  # Licensed under the Open Software License version 3.0
10      self.__matrix_type=0  # http://www.opensource.org/licenses/osl-3.0.php
11       def trace(self,text):  #
12           if self.__trace: print text  ########################################################
13    
14       def getDomain(self):  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15          return self.__domain  Earth Systems Science Computational Center (ESSCC)
16       def getTheta(self):  http://www.uq.edu.au/esscc
17          return self.__theta  Primary Business: Queensland, Australia"""
18       def getDt_max(self):  __license__="""Licensed under the Open Software License version 3.0
19          return self.__dt_max  http://www.opensource.org/licenses/osl-3.0.php"""
20       def getNumEquations(self):  __url__="https://launchpad.net/escript-finley"
21          return self.__num_equations  
22       def reduced(self):  from esys.escript import *
23           return False  from esys.escript.linearPDEs import LinearPDE, TransportPDE
      def getFunctionSpace(self):  
         if self.reduced():  
            return ReducedSolution(self.getDomain())  
         else:  
            return Solution(self.getDomain())  
   
   
      def __getNewTransportProblem(self):  
        """  
        returns an instance of a new operator  
        """  
        self.trace("New Transport problem is allocated.")  
        return self.getDomain().newTransportProblem( \  
                                self.getTheta(),  
                                self.getDt_max(),  
                                self.getNumEquations(), \  
                                self.getFunctionSpace(), \  
                                self.__matrix_type)  
   
      def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(),  
               d=Data(),y=Data(),d_contact=Data(),y_contact=Data()):  
          self.__transport_problem=self.__getNewTransportProblem()  
      if self.getNumEquations() ==1 :  
         self.__source=Data(0.0,(),self.getFunctionSpace())  
          else:  
              self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())  
      self.getDomain().addPDEToTransportProblem(  
                  self.__transport_problem,  
              self.__source,  
              M,A,B,C,D,X,Y,d,y,d_contact,y_contact)  
   
      def setInitialSolution(self,u):  
          self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))  
      def solve(self,dt):  
        return self.__transport_problem.solve(self.__source,dt,{"verbose" : True })  
24  from esys.finley import Rectangle  from esys.finley import Rectangle
25    
26  dom=Rectangle(40,20,l0=2)  # dom=Rectangle(12,8,l0=1.5)
27  fc=TransportPDE(dom,num_equations=1,theta=1.0,dt_max=2.5e-2/10)  # dom=Rectangle(24,16,l0=1.5)
28  fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])  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()  x=dom.getX()
37  x_0=[0.3,0.3]  x_0=[0.5,0.5]
38  sigma=0.08  sigma=0.075
39  u=1.  u0=1.
40  for i in range(dom.getDim()):  for i in range(dom.getDim()):
41      u=u*exp(-(x[i]-x_0[i])**2/sigma**2)      u0=u0*exp(-(x[i]-x_0[i])**2/sigma**2)
 u/=Lsup(u)  
42    
43  u=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)  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  c=0
51  saveVTK("u.%s.xml"%c,u=u)  saveVTK("u.%s.vtu"%c,u=u0)
52  fc.setInitialSolution(u)  fc.setInitialSolution(u0)
53    
54  dt=2.5e-2  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.  t=0.
61  while t<25*dt:  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_end and c< c_stop:
66      print "time step t=",t+dt        print "time step t=",t+dt  
67      u=fc.solve(dt)        u=fc.solve(dt)  
68      print "range u",inf(u),sup(u),integrate(u,Function(dom))      print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
69      c+=1      c+=1
     saveVTK("u.%s.xml"%c,u=u)  
70      t+=dt      t+=dt
71        if 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_end and c< c_stop:
87           print "time step t=",t+dt    
88           pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
89           u=pde.getSolution(verbose=True)
90           print t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom))
91           c+=1
92           t+=dt
93           if 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)

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

  ViewVC Help
Powered by ViewVC 1.1.26