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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1401 - (show annotations)
Fri Jan 25 04:31:18 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 5363 byte(s)
rewrite antidiffusion calculation to avoid coloring for OPENMP parallelization
1 #
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 from esys.escript import *
55 from esys.escript.linearPDEs import LinearSinglePDE
56
57 class TransportPDE(object):
58 def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
59 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 self.__trace=trace
65 self.__matrix_type=0
66 def trace(self,text):
67 if self.__trace: print text
68
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 return self.__num_equations
77 def reduced(self):
78 return False
79 def getFunctionSpace(self):
80 if self.reduced():
81 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 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 self.__transport_problem=self.__getNewTransportProblem()
101 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
110 def setInitialSolution(self,u):
111 self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
112 def solve(self,dt):
113 return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : 1.e-8})
114
115
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 x=dom.getX()
127 u0=1/(4.*pi*E*T0)**(DIM/2.)*exp(-length(dom.getX()-getCenter(T0))**2/(4.*E*T0))
128
129 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 c=0
149 saveVTK("u.%s.xml"%c,u=u0)
150 fc.setInitialSolution(u0)
151 dt=1.266539e-02*10
152 t=T0
153 while t<T0+2:
154 print "time step t=",t+dt
155 u=fc.solve(dt)
156 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 c+=1
168 t+=dt
169 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