/[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 1400 - (hide annotations)
Thu Jan 24 06:04:31 2008 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 5364 byte(s)
better test example for upwinding added
1 gross 1400 #
2     # upwinding test moving a Gaussian hill around
3     #
4     # we solve U_,t - E *u_,ii + v_i u_,i =0 (E is small)
5     #
6     # the solution is given as u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) )
7     #
8     # where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and
9     #
10     # x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3
11     #
12     # the solution is started from some time T0>0.
13     #
14     # We are using five quality messurements for u_h
15     #
16     # - inf(u_h) > 0
17     # - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1
18     # - integrate(u_h) ~ 1
19     # - | x_0h-x_0 | ~ 0 where x_0h = integrate(x*u_h)
20     # - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )
21     #
22     #
23     from math import pi, ceil
24     NE=128/2
25     DIM=2
26     THETA=0.5
27     OMEGA0=1.
28     ALPHA=pi/4
29     T0=0.5*pi
30     E=5e-4 * 2
31     TEST_SUPG=False
32    
33     def getCenter(t):
34     if DIM==2:
35     x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5]
36     else:
37     x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5]
38     return x0
39     def QUALITY(t,u_h):
40     dom=u_h.getDomain()
41     x=dom.getX()
42     a=inf(u_h)
43     b=sup(u_h)*(4*pi*E*t)**(DIM/2.)-1.
44     c=integrate(u_h,Function(dom))-1.
45     x0=getCenter(t)
46     x0h=integrate(x*u_h,Function(dom))
47     d=length(x0-x0h)
48     sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom)))
49     if DIM == 3: sigma_h2*=sqrt(2./3.)
50     e=sigma_h2/sqrt(4*E*t)-1
51     return a,b,c,d,e
52    
53    
54 gross 1368 from esys.escript import *
55 gross 1400 from esys.escript.linearPDEs import LinearSinglePDE
56 gross 1368
57     class TransportPDE(object):
58 gross 1370 def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
59 gross 1368 self.__domain=domain
60     self.__num_equations=num_equations
61     self.__theta=theta
62     self.__dt_max=dt_max
63     self.__transport_problem=None
64 gross 1370 self.__trace=trace
65     self.__matrix_type=0
66     def trace(self,text):
67     if self.__trace: print text
68 gross 1368
69     def getDomain(self):
70     return self.__domain
71     def getTheta(self):
72     return self.__theta
73     def getDt_max(self):
74     return self.__dt_max
75     def getNumEquations(self):
76 gross 1370 return self.__num_equations
77     def reduced(self):
78     return False
79 gross 1368 def getFunctionSpace(self):
80 gross 1370 if self.reduced():
81 gross 1368 return ReducedSolution(self.getDomain())
82     else:
83     return Solution(self.getDomain())
84    
85    
86     def __getNewTransportProblem(self):
87     """
88     returns an instance of a new operator
89     """
90     self.trace("New Transport problem is allocated.")
91     return self.getDomain().newTransportProblem( \
92     self.getTheta(),
93     self.getDt_max(),
94     self.getNumEquations(), \
95     self.getFunctionSpace(), \
96     self.__matrix_type)
97    
98 gross 1370 def setValue(self,M=Data(),A=Data(),B=Data(),C=Data(),D=Data(),X=Data(),Y=Data(),
99     d=Data(),y=Data(),d_contact=Data(),y_contact=Data()):
100 gross 1368 self.__transport_problem=self.__getNewTransportProblem()
101 gross 1370 if self.getNumEquations() ==1 :
102     self.__source=Data(0.0,(),self.getFunctionSpace())
103     else:
104     self.__source=Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
105     self.getDomain().addPDEToTransportProblem(
106     self.__transport_problem,
107     self.__source,
108     M,A,B,C,D,X,Y,d,y,d_contact,y_contact)
109 gross 1368
110 gross 1370 def setInitialSolution(self,u):
111     self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
112     def solve(self,dt):
113 gross 1400 return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : 1.e-8})
114 gross 1370
115 gross 1400
116     from esys.finley import Rectangle, Brick
117    
118    
119     if DIM==2:
120     dom=Rectangle(NE,NE)
121     else:
122     dom=Brick(NE,NE,NE)
123     dom.setX(2*dom.getX()-1)
124    
125     # set initial value
126 gross 1370 x=dom.getX()
127 gross 1400 u0=1/(4.*pi*E*T0)**(DIM/2.)*exp(-length(dom.getX()-getCenter(T0))**2/(4.*E*T0))
128 gross 1371
129 gross 1400 print "QUALITY ",QUALITY(T0,u0)
130    
131     x=Function(dom).getX()
132     if DIM == 2:
133     V=OMEGA0*(x[0]*[0,1]+x[1]*[-1,0])
134     else:
135     V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
136     #===================
137     fc=TransportPDE(dom,num_equations=1,theta=THETA)
138     x=Function(dom).getX()
139     fc.setValue(M=Scalar(1.,Function(dom)),C=V,A=-Scalar(E,Function(dom))*kronecker(dom))
140     #==============
141     if TEST_SUPG:
142     supg=LinearSinglePDE(dom)
143     supg.setValue(D=1.)
144     supg.setSolverMethod(supg.LUMPING)
145     dt_supg=1./(1./inf(dom.getSize()/length(V))+1./inf(dom.getSize()**2/E))*0.3
146     u_supg=u0*1.
147    
148 gross 1370 c=0
149 gross 1400 saveVTK("u.%s.xml"%c,u=u0)
150     fc.setInitialSolution(u0)
151     dt=1.266539e-02*10
152     t=T0
153     while t<T0+10:
154 gross 1370 print "time step t=",t+dt
155     u=fc.solve(dt)
156 gross 1400 if TEST_SUPG:
157     #========== supg tests ================
158     nn=max(ceil(dt/dt_supg),1.)
159     dt2=dt/nn
160     nnn=0
161     while nnn<nn :
162     supg.setValue(X=-dt2/2*E*grad(u_supg),Y=u_supg+dt2/2*inner(V,grad(u_supg)))
163     u2=supg.getSolution()
164     supg.setValue(X=-dt2*E*grad(u2),Y=u_supg+dt2*inner(V,grad(u2)))
165     u_supg=supg.getSolution()
166     nnn+=1
167 gross 1370 c+=1
168     t+=dt
169 gross 1400 print "QUALITY FCT: ",QUALITY(t,u)
170     if TEST_SUPG:
171     print "QUALITY SUPG: ",QUALITY(t,u_supg)
172     # saveVTK("u.%s.xml"%c,u=u,u_supg=u_supg)
173     else:
174     # saveVTK("u.%s.xml"%c,u=u)
175     pass
176     # if c == 20: 1/0

  ViewVC Help
Powered by ViewVC 1.1.26