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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3033 - (show annotations)
Wed Jun 2 10:17:00 2010 UTC (9 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 8487 byte(s)
fix for overflow
1 # -*- coding: utf-8 -*-
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 #
23 # Flux corrected transport solver benchmark
24 # we are moving a Gaussian hill around
25 #
26 # we solve a* U_{,t} - b *u_{,ii} + c_i u_{,i} + (d_i* u)_{,i}}=0
27 #
28 # U(0)= U0 * exp ( - |x-x_0(t)|^2/(4*s) )
29 #
30 # with a>0, b>=0, s>0
31 #
32 # we set E=b/a v=c/a w=d/a
33 #
34 # the solution is given as u(x,t)=U0*s^{dim/2}/(s+E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*(s+E*t)) )
35 #
36 # with x_0(t) = X0 + (v+w)*t
37 #
38 #
39 # |x-x_0(t)|^2/(4*(s+E*t)) < - log(TAU) for all time and all x in the
40 # domain where TAU is a given relative tolerance
41 #
42 #
43 # this holds if
44 #
45 # |x_i-X0_i-(v_i+w_i)*tend | < - log(TAU) * 4*(s+E*tend)
46 # |x_i-X0_i | < - log(TAU) * 4*s
47 #
48 # X0_i=min(- log(TAU) * 4*s,-(v_i+w_i)*tend - log(TAU) * 4*(s+E*tend))
49 # l_i=X0_i+min((v_i+w_i)*tend - log(TAU) * 4*(s+E*tend)) , - log(TAU) * 4*s)
50 #
51 DIM=2
52 NE_MAX=35000
53 VERBOSITY=True
54 TOL=1.e-8
55 TAU=1e-10
56 VTK_DIR="output"
57 S_RAISE=1.5
58 S_MIN=0
59 S_MAX=0
60 TABS = [ 'dx', 'dt', 'peclet', 'error', 'sup', 'integral', 'center', 'width', 'time' ]
61 from math import pi, ceil
62 from time import time as clock
63 from esys.finley import Rectangle, Brick
64 from esys.escript import *
65 from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
66 S_MAX=-0.5/log(TAU)/4
67
68
69 mkDir(VTK_DIR)
70
71 def uRef(dom,t,E,s,v,x0, onElements=False):
72 if onElements:
73 x=Function(dom).getX()
74 else:
75 x=dom.getX()
76 X=x0[:dom.getDim()]+v[:dom.getDim()]*t
77 u=(s/(s+E*t))**(dom.getDim()/2.) * exp(-length(x-X)**2/(4*(s+E*t)))
78 return u
79
80 NE_test = [ 10, 20, 40, 80, 160 ]
81 NE_test = [ 10, 20, 40, 80, 160, 320]
82 NE_test = [ 320 *2 ]
83 b_test=[ 0., 0.001, 0.1, 1., 10, 100.]
84 b_test=[ 1.0 *0]
85 c_test=[ 0., 0.001, 0.1, 1., 10, 100. ]
86 c_test=[ 100.]
87 d_test=[ 0.0 ]
88 s_test=[0.0005 ]
89 # s=[ 0.001, 0.003, 0.01 ]
90
91 a=1.
92 X0=[0.5,0.5,0.5]
93 step=0.1
94 s_max=0.01 # don't use an s which is bigger than this as the profile will not fit into the domian
95
96
97 def getDirection(dim, d="x"):
98 k=kronecker(dim)
99 if d=="x":
100 return k[0]
101 elif d=="y":
102 return k[1]
103 elif d=="z" and dim>2:
104 return k[2]
105 elif d=="xy":
106 return (k[0]+k[1])/sqrt(2.)
107 elif d=="yz":
108 return (k[1]+k[2])/sqrt(2.)
109 elif d=="zx" and dim>2:
110 return (k[2]+k[0])/sqrt(2.)
111 elif d=="xyz" and dim>2:
112 return (k[1]+k[2]+k[0])/sqrt(3.)
113 else:
114 raise ValueError,"Cannot identify direction %s"%d
115
116 #================
117 def QUALITY(u_h,u_ref):
118 dom=u_h.getDomain()
119 x=dom.getX()
120 out = {}
121 out["sup"]=abs(sup(u_h)-sup(u_ref))/abs(sup(u_ref))
122 i_h=integrate(u_h)
123 i_ref=integrate(u_ref)
124 out["integral"]=abs(i_h-i_ref)/abs(i_ref)
125 x0h=integrate(x*u_h)/i_h
126 x0ref=integrate(x*u_ref)/i_ref
127 out["center"]=length(x0h-x0ref)/(1.+step)
128
129 sigma_h=sqrt(integrate(length(x-x0h)**2*u_h))
130 sigma_ref=sqrt(integrate(length(x-x0ref)**2*u_ref))
131 out["width"]=abs(sigma_h-sigma_ref)/abs(sigma_ref)
132 out["error"]=sqrt(integrate((u_ref-u_h)**2))
133
134 return out
135
136
137
138 def XXX(dim, s, h,b,c,d,c_dir="x", d_dir="x", a=1.):
139 v_c=c/a*getDirection(dim,c_dir)
140 v_d=d/a*getDirection(dim,d_dir)
141 v = (v_c+v_d)
142 E=b/a
143 if abs(E)>0:
144 if length(v)>0:
145 tend=min(s/E, s/length(v))
146 else:
147 tend=s/E
148 else:
149 if length(v)>0:
150 tend=s/length(v)
151 else:
152 tend=0
153 tend*=(1.-S_RAISE)
154 if VERBOSITY:
155 print "="*100
156 print "Start test dim = %d , h=%e, b=%e, c=%e, d=%e, c_dir=%s, d_dir=%s, a=%e, s=%e"%(dim, h,b,c,d,c_dir, d_dir, a, s)
157 print "="*100
158 print "initial width s = ",s
159 print "diffusion = ",E
160 print "total velocity = ",v
161 print "tend = ", tend
162 print "tolerance = ",TOL
163 print "XX number of elements over s =",s/h
164 X0_0=min(- log(TAU) * 4*s,-v[0]*tend - log(TAU) * 4*(s+E*tend))
165 X0_1=min(- log(TAU) * 4*s,-v[1]*tend - log(TAU) * 4*(s+E*tend))
166 l_0=X0_0+min(v[0]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)
167 l_1=X0_1+min(v[1]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)
168 NE_0=max(int(l_0/h+0.5),1)
169 NE_1=max(int(l_1/h+0.5),1)
170 if dim == 2:
171 if VERBOSITY: print "%d x %d grid over %e x %e with element size %e."%(NE_0,NE_1,l_0,l_1,h)
172 if NE_0*NE_1 > NE_MAX:
173 raise ValueError,"too many elements %s."%(NE_0*NE_1,)
174 dom=Rectangle(n0=NE_0,n1=NE_1,l0=l_0,l1=l_1)
175 x0=[X0_0, X0_1]
176 else:
177 X0_2=min(- log(TAU) * 4*s,-v[2]*tend - log(TAU) * 4*(s+E*tend))
178 l_2=X0_2+min(v[2]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)
179 NE_2=max(int(l_2/h+0.5),1)
180 if VERBOSITY: print "%d x %d x %d grid over %e x %e x %e with element size %e."%(NE_0,NE_1,NE_2,l_0,l_1,l_2,h)
181 if NE_0*NE_1*NE_2 > NE_MAX:
182 raise ValueError,"too many elements %s."%(NE_0*NE_1*NE_2,)
183 dom=Brick(n0=NE_0,n1=NE_1, ne2=NE_2, l0=l_0,l1=l_1, l2=l_2)
184 x0=[X0_0, X0_1, X0_2]
185 if VERBOSITY:
186 print "initial location = ",x0
187
188 fc_BE=TransportPDE(dom,numEquations=1,useBackwardEuler=True)
189 fc_BE.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
190 fc_BE.getSolverOptions().setVerbosity(VERBOSITY)
191 fc_BE.getSolverOptions().setTolerance(TOL)
192 #
193 # fc_BE.getSolverOptions().setPreconditioner(fc_BE.getSolverOptions().GAUSS_SEIDEL)
194 fc_BE.getSolverOptions().setNumSweeps(1)
195 if VERBOSITY: print "Backward Euler Transport created"
196
197 fc_CN=TransportPDE(dom,numEquations=1,useBackwardEuler=False)
198 fc_CN.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
199 fc_CN.getSolverOptions().setVerbosity(VERBOSITY)
200 fc_CN.getSolverOptions().setTolerance(TOL)
201
202 #fc_CN.getSolverOptions().setPreconditioner(fc_CN.getSolverOptions().GAUSS_SEIDEL)
203 fc_CN.getSolverOptions().setNumSweeps(2)
204 if VERBOSITY: print "Crank Nicolson Transport created"
205 dt_CN=fc_CN.getSafeTimeStepSize()
206 if VERBOSITY: print "time step size by Crank Nicolson=",dt_CN
207 if tend>0:
208 dt=main(dt_CN, tend)
209 else:
210 dt=1.
211 U0=uRef(dom,0,E,s,v,X0)
212 U0_e=uRef(dom,0,E,s,v,X0,True)
213 fc_CN.setInitialSolution(U0)
214 fc_BE.setInitialSolution(U0)
215 initial_error_L2=sqrt(integrate((U0-U0_e)**2))
216 if VERBOSITY:
217 print "used time step size =",dt
218 print "initial integral = ",integrate(U0_e)
219 print "initial error = ",initial_error_L2
220
221
222
223 return
224 if False:
225 # dt2=min(dt2,20*dt)
226 if VERBOSITY: print "Solve Backward Euler:"
227 t0=clock()
228 u_BE=fc_BE.getSolution(10.*dt)
229 u_BE=fc_BE.getSolution(10.*dt)
230 t0=clock()-t0
231 out=QUALITY(u_BE,uRef(dom,dt2,E,s,v,X0))
232 #
233 # print integrate(U0), U0
234 # print integrate(u_BE), u_BE
235 # print integrate(uRef(dom,dt2,E,s,v,X0)), uRef(dom,dt2,E,s,v,X0)
236 # saveVTK("bb_be.vtu",u0=U0,u_BE=u_BE, uRef=uRef(dom,dt2,E,s,v,X0) )
237 # 1/0
238 #
239 else:
240 if VERBOSITY: print "Solve Crank Nicolson:"
241 t0=clock()
242 u_CN=fc_CN.getSolution(dt)
243 t0=clock()-t0
244 out=QUALITY(u_CN,uRef(dom,dt,E,s,v,X0))
245 out['time']=t0/NN
246 out['dt']=dt
247 out['dx']=1./NE
248 if abs(b)>0:
249 out["peclet"]=length(v)/(b*NE)
250 else:
251 out["peclet"]=9999999.
252 # saveVTK("bb.vtu",u0=U0,u_CN=u_CN, uRef=uRef(dom,dt2,E,s,v,X0) )
253 return out
254
255
256 for s_fac in [ 1., 0.5, 0.25, 0.125 ] :
257 for h_fac in [1., 0.5, 0.25, 0.125 ]:
258 XXX(DIM,s=S_MAX*s_fac,h=10*S_MAX*s_fac*h_fac,b=0,c=0,d=0,c_dir="x", d_dir="x")
259 1/0
260
261 txt=""
262 for t in TABS: txt+=" "+t
263 for b in b_test:
264 for c in c_test:
265 for d in d_test:
266 for s in s_test:
267 for NE in NE_test:
268 dat=XXX(NE, a,b,c,d,c_dir="x", d_dir="x",s=s)
269 txt+="\n"
270 for t in TABS: txt+=" %3.2e"%dat[t]
271
272 print txt

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26