/[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 1407 - (show annotations)
Mon Feb 4 06:45:48 2008 UTC (11 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 2776 byte(s)
new upwinding algorithm (still fails)
1 from esys.escript import *
2
3 class TransportPDE(object):
4 def __init__(self,domain,num_equations=1,theta=0.,trace=True):
5 self.__domain=domain
6 self.__num_equations=num_equations
7 self.__theta=theta
8 self.__transport_problem=None
9 self.__trace=trace
10 self.__matrix_type=0
11 def trace(self,text):
12 if self.__trace: print text
13
14 def getDomain(self):
15 return self.__domain
16 def getTheta(self):
17 return self.__theta
18 def getNumEquations(self):
19 return self.__num_equations
20 def reduced(self):
21 return False
22 def getFunctionSpace(self):
23 if self.reduced():
24 return ReducedSolution(self.getDomain())
25 else:
26 return Solution(self.getDomain())
27
28 def getSafeTimeStepSize(self):
29 return self.__transport_problem.getSafeTimeStepSize()
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.getNumEquations(), \
40 self.getFunctionSpace(), \
41 self.__matrix_type)
42
43 def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(),
44 d=Data(),y=Data(),d_contact=Data(),y_contact=Data()):
45 self.__transport_problem=self.__getNewTransportProblem()
46 if self.getNumEquations() ==1 :
47 self.__source=Data(0.0,(),self.getFunctionSpace())
48 else:
49 self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
50 self.getDomain().addPDEToTransportProblem(
51 self.__transport_problem,
52 self.__source,
53 M,A,B,C,D,X,Y,d,y,d_contact,y_contact)
54
55 def setInitialSolution(self,u):
56 self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
57 def solve(self,dt):
58 return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : 1.e-6})
59 from esys.finley import Rectangle
60
61 dom=Rectangle(6,3,l0=1.5)
62 # dom=Rectangle(120,80,l0=1.5)
63 fc=TransportPDE(dom,num_equations=1,theta=0.5)
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=u*exp(-(x[i]-x_0[i])**2/sigma**2)
71 u/=Lsup(u)
72
73 u=whereNonPositive(abs(x[0]-0.4)-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<50*dt:
81 print "time step t=",t+dt
82 u=fc.solve(dt)
83 print "range u",inf(u),sup(u),integrate(u,Function(dom))
84 c+=1
85 saveVTK("u.%s.xml"%c,u=u)
86 t+=dt
87 if c == 20: 1/0

  ViewVC Help
Powered by ViewVC 1.1.26