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

Annotation of /trunk/finley/test/python/FCT_test2.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1388 - (hide annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 8 months ago) by trankine
Original Path: trunk/finley/test/python/tp.py
File MIME type: text/x-python
File size: 2739 byte(s)
And get the *(&(*&(* name right
1 gross 1368 from esys.escript import *
2    
3     class TransportPDE(object):
4 gross 1370 def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
5 gross 1368 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 gross 1370 self.__trace=trace
11     self.__matrix_type=0
12     def trace(self,text):
13     if self.__trace: print text
14 gross 1368
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 gross 1370 return self.__num_equations
23     def reduced(self):
24     return False
25 gross 1368 def getFunctionSpace(self):
26 gross 1370 if self.reduced():
27 gross 1368 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 gross 1370 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 gross 1368 self.__transport_problem=self.__getNewTransportProblem()
47 gross 1370 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 gross 1368
56 gross 1370 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 gross 1371 dom=Rectangle(20,20)
63 gross 1370 fc=TransportPDE(dom,num_equations=1,theta=0.0,dt_max=0.)
64 gross 1371 fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])
65 gross 1370 x=dom.getX()
66 gross 1371 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 gross 1375 # u=whereNonPositive(abs(x[0]-0.3)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)
74 gross 1370 c=0
75     saveVTK("u.%s.xml"%c,u=u)
76     fc.setInitialSolution(u)
77 gross 1371
78     dt=2.5e-2
79 gross 1370 t=0.
80 gross 1371 while t<15*dt:
81 gross 1370 print "time step t=",t+dt
82     u=fc.solve(dt)
83 gross 1375 print "range u",inf(u),sup(u),integrate(u,Function(dom))
84 gross 1370 c+=1
85     saveVTK("u.%s.xml"%c,u=u)
86     t+=dt

  ViewVC Help
Powered by ViewVC 1.1.26