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

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