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

revision 1392 by gross, Mon Jan 21 06:20:54 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 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-2009 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
11       def trace(self,text):  #
12           if self.__trace: print text  ########################################################
13
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"""
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.__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