/[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 1407 by gross, Mon Feb 4 06:45:48 2008 UTC revision 1410 by gross, Thu Feb 7 04:24:00 2008 UTC
# Line 1  Line 1 
1  from esys.escript import *  from esys.escript import *
2    from esys.escript.linearPDEs import LinearPDE, TransportPDE
 class TransportPDE(object):  
      def __init__(self,domain,num_equations=1,theta=0.,trace=True):  
         self.__domain=domain  
         self.__num_equations=num_equations  
         self.__theta=theta  
         self.__transport_problem=None  
     self.__trace=trace  
     self.__matrix_type=0  
      def trace(self,text):  
          if self.__trace: print text  
   
      def getDomain(self):  
         return self.__domain  
      def getTheta(self):  
         return self.__theta  
      def getNumEquations(self):  
         return self.__num_equations  
      def reduced(self):  
          return False  
      def getFunctionSpace(self):  
         if self.reduced():  
            return ReducedSolution(self.getDomain())  
         else:  
            return Solution(self.getDomain())  
   
      def getSafeTimeStepSize(self):  
         return self.__transport_problem.getSafeTimeStepSize()  
   
   
      def __getNewTransportProblem(self):  
        """  
        returns an instance of a new operator  
        """  
        self.trace("New Transport problem is allocated.")  
        return self.getDomain().newTransportProblem( \  
                                self.getTheta(),  
                                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 , "tolerance" : 1.e-6})  
3  from esys.finley import Rectangle  from esys.finley import Rectangle
4    
5  dom=Rectangle(6,3,l0=1.5)  # 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  # dom=Rectangle(120,80,l0=1.5)  # dom=Rectangle(120,80,l0=1.5)
10  fc=TransportPDE(dom,num_equations=1,theta=0.5)  V=Scalar(1.,Function(dom))*[-1.,0]
11  fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])  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  x=dom.getX()  x=dom.getX()
16  x_0=[0.3,0.3]  x_0=[0.5,0.5]
17  sigma=0.08  sigma=0.075
18  u=1.  u0=1.
19  for i in range(dom.getDim()):  for i in range(dom.getDim()):
20      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)  
21    
22  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)
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  c=0  c=0
30  saveVTK("u.%s.xml"%c,u=u)  saveVTK("u.%s.xml"%c,u=u0)
31  fc.setInitialSolution(u)  fc.setInitialSolution(u0)
32    
33  dt=2.5e-2  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  t=0.  t=0.
40  while t<50*dt:  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_end and c< c_stop:
45      print "time step t=",t+dt        print "time step t=",t+dt  
46      u=fc.solve(dt)        u=fc.solve(dt)  
47      print "range u",inf(u),sup(u),integrate(u,Function(dom))      print t+dt,": range u",inf(u),sup(u),integrate(u,Function(dom))
48      c+=1      c+=1
     saveVTK("u.%s.xml"%c,u=u)  
49      t+=dt      t+=dt
50      if c == 20: 1/0      if 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_end and c< c_stop:
66           print "time step t=",t+dt    
67           pde.setValue(Y=u+(1.-THETA)*dt*inner(V,grad(u)))
68           u=pde.getSolution(verbose=True)
69           print t+dt,": range u2",inf(u),sup(u),integrate(u,Function(dom))
70           c+=1
71           t+=dt
72           if 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)
Legend:
Removed from v.1407  
changed lines
  Added in v.1410
  ViewVC Help
Powered by ViewVC 1.1.26