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

Contents of /trunk/finley/test/python/tp.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1370 - (show annotations)
Wed Jan 2 09:21:43 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 2597 byte(s)
explicit upwinding scheme added.
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(10,10)
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 u=whereNonPositive(abs(x[0]-0.3)-0.1)*whereNonPositive(abs(x[1]-0.3)-0.1)
67 c=0
68 saveVTK("u.%s.xml"%c,u=u)
69 fc.setInitialSolution(u)
70 dt=0.301
71 t=0.
72 while t<dt:
73 print "time step t=",t+dt
74 u=fc.solve(dt)
75 print "range u",inf(u),sup(u)
76 c+=1
77 saveVTK("u.%s.xml"%c,u=u)
78 t+=dt

  ViewVC Help
Powered by ViewVC 1.1.26