test/python/tp.py

from esys.escript import *

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,{})
from esys.finley import Rectangle

dom=Rectangle(20,20)
fc=TransportPDE(dom,num_equations=1,theta=0.0,dt_max=0.)
fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])
x=dom.getX()
x_0=[0.3,0.3]
sigma=0.08
u=1.
for i in range(dom.getDim()):
    u=u*exp(-(x[i]-x_0[i])**2/sigma**2)
u/=Lsup(u)

u=whereNonPositive(abs(x[0]-0.3)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
c=0
saveVTK("u.%s.xml"%c,u=u)
fc.setInitialSolution(u)

dt=2.5e-2
t=0.
while t<15*dt:
    print "time step t=",t+dt   
    u=fc.solve(dt)  
    print "range u",inf(u),sup(u)
    c+=1
    saveVTK("u.%s.xml"%c,u=u)
    t+=dt
1	from esys.escript import *
2
3	class TransportPDE(object):
4	def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
5	self.__domain=domain
6	self.__num_equations=num_equations
7	self.__theta=theta
8	self.__dt_max=dt_max
9	self.__transport_problem=None
10	self.__trace=trace
11	self.__matrix_type=0
12	def trace(self,text):
13	if self.__trace: print text
14
15	def getDomain(self):
16	return self.__domain
17	def getTheta(self):
18	return self.__theta
19	def getDt_max(self):
20	return self.__dt_max
21	def getNumEquations(self):
22	return self.__num_equations
23	def reduced(self):
24	return False
25	def getFunctionSpace(self):
26	if self.reduced():
27	return ReducedSolution(self.getDomain())
28	else:
29	return Solution(self.getDomain())
30
31
32	def __getNewTransportProblem(self):
33	"""
34	returns an instance of a new operator
35	"""
36	self.trace("New Transport problem is allocated.")
37	return self.getDomain().newTransportProblem( \
38	self.getTheta(),
39	self.getDt_max(),
40	self.getNumEquations(), \
41	self.getFunctionSpace(), \
42	self.__matrix_type)
43
44	def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(),
45	d=Data(),y=Data(),d_contact=Data(),y_contact=Data()):
46	self.__transport_problem=self.__getNewTransportProblem()
47	if self.getNumEquations() ==1 :
48	self.__source=Data(0.0,(),self.getFunctionSpace())
49	else:
50	self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
51	self.getDomain().addPDEToTransportProblem(
52	self.__transport_problem,
53	self.__source,
54	M,A,B,C,D,X,Y,d,y,d_contact,y_contact)
55
56	def setInitialSolution(self,u):
57	self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
58	def solve(self,dt):
59	return self.__transport_problem.solve(self.__source,dt,{})
60	from esys.finley import Rectangle
61
62	dom=Rectangle(20,20)
63	fc=TransportPDE(dom,num_equations=1,theta=0.0,dt_max=0.)
64	fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])
65	x=dom.getX()
66	x_0=[0.3,0.3]
67	sigma=0.08
68	u=1.
69	for i in range(dom.getDim()):
70	u=uexp(-(x[i]-x_0[i])2/sigma*2)
71	u/=Lsup(u)
72
73	u=whereNonPositive(abs(x[0]-0.3)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
74	c=0
75	saveVTK("u.%s.xml"%c,u=u)
76	fc.setInitialSolution(u)
77
78	dt=2.5e-2
79	t=0.
80	while t<15*dt:
81	print "time step t=",t+dt
82	u=fc.solve(dt)
83	print "range u",inf(u),sup(u)
84	c+=1
85	saveVTK("u.%s.xml"%c,u=u)
86	t+=dt