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

revision 1409 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.__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