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

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

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

revision 1401 by gross, Fri Jan 25 04:31:18 2008 UTC revision 1410 by gross, Thu Feb 7 04:24:00 2008 UTC
# Line 20  Line 20 
20  #     - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )  #     - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )
21  #  #
22  #  #
23    from esys.escript import *
24    from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
25    from esys.finley import Rectangle, Brick
26  from math import pi, ceil  from math import pi, ceil
27  NE=128/2  NE=128
28  DIM=2  DIM=2
29  THETA=0.5  THETA=0.5
30  OMEGA0=1.  OMEGA0=1.
31  ALPHA=pi/4  ALPHA=pi/4
32  T0=0.5*pi  T0=0.5*pi
33  E=5e-4 * 2  T_END=2.5*pi
34  TEST_SUPG=False  dt=1e-3*10
35    E=1.e-3
36    TEST_SUPG=False # or True
37    
38    
39  def getCenter(t):  def getCenter(t):
40     if DIM==2:     if DIM==2:
41        x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5]        x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5]
42          x0=[-sin(OMEGA0*t)*0.5,cos(OMEGA0*t)*0.5]
43     else:     else:
44        x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5]        x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5]
45     return x0     return x0
# Line 48  def QUALITY(t,u_h): Line 55  def QUALITY(t,u_h):
55       sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom)))       sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom)))
56       if DIM == 3: sigma_h2*=sqrt(2./3.)       if DIM == 3: sigma_h2*=sqrt(2./3.)
57       e=sigma_h2/sqrt(4*E*t)-1                   e=sigma_h2/sqrt(4*E*t)-1            
58       return a,b,c,d,e       # return a,b,c,e,1./(4*pi*E*t)**(DIM/2.)
59         return b, d
60         # return a,b,c,d,e
61                
62    
 from esys.escript import *  
 from esys.escript.linearPDEs import LinearSinglePDE  
   
 class TransportPDE(object):  
      def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):  
         self.__domain=domain  
         self.__num_equations=num_equations  
         self.__theta=theta  
         self.__dt_max=dt_max  
         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 getDt_max(self):  
         return self.__dt_max  
      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 __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 , "tolerance" : 1.e-8})  
   
   
 from esys.finley import Rectangle, Brick  
63    
64    
65  if DIM==2:  if DIM==2:
# Line 130  print "QUALITY ",QUALITY(T0,u0) Line 76  print "QUALITY ",QUALITY(T0,u0)
76    
77  x=Function(dom).getX()  x=Function(dom).getX()
78  if DIM == 2:  if DIM == 2:
79     V=OMEGA0*(x[0]*[0,1]+x[1]*[-1,0])     V=OMEGA0*(x[0]*[0,-1]+x[1]*[1,0])
80  else:  else:
81     V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])     V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
82  #===================  #===================
# Line 148  if TEST_SUPG: Line 94  if TEST_SUPG:
94  c=0  c=0
95  saveVTK("u.%s.xml"%c,u=u0)  saveVTK("u.%s.xml"%c,u=u0)
96  fc.setInitialSolution(u0)  fc.setInitialSolution(u0)
 dt=1.266539e-02*10  
97  t=T0  t=T0
98  while t<T0+2:  while t<T_END:
99      print "time step t=",t+dt        print "time step t=",t+dt  
100      u=fc.solve(dt)        u=fc.solve(dt)  
101      if TEST_SUPG:      if TEST_SUPG:
# Line 166  while t<T0+2: Line 111  while t<T0+2:
111              nnn+=1              nnn+=1
112      c+=1      c+=1
113      t+=dt      t+=dt
114      print "QUALITY FCT: ",QUALITY(t,u)      print "QUALITY FCT: time = %s pi"%(t/pi),QUALITY(t,u)
115      if TEST_SUPG:      if TEST_SUPG:
116         print "QUALITY SUPG: ",QUALITY(t,u_supg)         print "QUALITY SUPG: ",QUALITY(t,u_supg)
117         # saveVTK("u.%s.xml"%c,u=u,u_supg=u_supg)         # saveVTK("u.%s.xml"%c,u=u,u_supg=u_supg)

Legend:
Removed from v.1401  
changed lines
  Added in v.1410

  ViewVC Help
Powered by ViewVC 1.1.26