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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 8843 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1 # -*- coding: utf-8 -*-
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 #
27 # Flux corrected transport solver benchmark
28 # we are moving a Gaussian hill around
29 #
30 # we solve a* U_{,t} - b *u_{,ii} + c_i u_{,i} + (d_i* u)_{,i}}=0
31 #
32 # U(0)= U0 * exp ( - |x-x_0(t)|^2/(4*s**2) )
33 #
34 # with a>0, b>=0, s>0
35 #
36 # we set E=b/a v=c/a w=d/a
37 #
38 # the solution is given as u(x,t)=U0*s^dim/(s**2+E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*(s**2+E*t)) )
39 #
40 # with x_0(t) = X0 + (v+w)*t
41 #
42 #
43 # the region |x-x_0(t)|^2/(4*(s**2+E*t)) < - log(TAU) is within the domain for all time
44 #
45 #
46 # this holds if
47 #
48 # |x_i-X0_i-(v_i+w_i)*tend | < sqrt(- log(TAU) * 4*(s**2+E*tend))=b0 and
49 # |x_i-X0_i | < sqrt(- log(TAU)) * 2*s = b1 implies 0<=x_i<=l_i
50 #
51 from math import pi, ceil
52 from time import time as clock
53 from esys.dudley import Rectangle, Brick
54 from esys.escript import *
55 from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
56 #
57 DIM=2
58 NE_MAX=300000
59 VERBOSITY=True
60 TOL=1.e-8
61 TAU=1e-10
62 VTK_DIR="output"
63
64 #==================
65 S_MIN=0
66 TABS = [ 'dx', 'dt', 'peclet', 'error', 'sup', 'integral', 'center', 'width', 'time' ]
67
68 a=1.
69 #==================
70
71
72 mkDir(VTK_DIR)
73
74 def uRef(dom,t,E,s,v,x0, onElements=False):
75 if onElements:
76 x=Function(dom).getX()
77 else:
78 x=dom.getX()
79 X=x0[:dom.getDim()]+v[:dom.getDim()]*t
80 u=(s**2/(s**2+E*t))**(dom.getDim()/2.) * exp(-length(x-X)**2/(4*(s**2+E*t)))
81 return u
82
83
84 def getDirection(dim, d="x"):
85 k=kronecker(dim)
86 if d=="x":
87 return k[0]
88 elif d=="y":
89 return k[1]
90 elif d=="z" and dim>2:
91 return k[2]
92 elif d=="xy":
93 return (k[0]+k[1])/sqrt(2.)
94 elif d=="yz":
95 return (k[1]+k[2])/sqrt(2.)
96 elif d=="zx" and dim>2:
97 return (k[2]+k[0])/sqrt(2.)
98 elif d=="xyz" and dim>2:
99 return (k[1]+k[2]+k[0])/sqrt(3.)
100 else:
101 raise ValueError("Cannot identify direction %s"%d)
102
103 def QUALITY(u_h,u_ref):
104 u_h_e=interpolate(u_h,u_ref.getFunctionSpace())
105 x=u_ref.getFunctionSpace().getX()
106 out = {}
107 out["error"]=sqrt(integrate((u_ref-u_h_e)**2))/sqrt(integrate(u_ref**2))
108 out["sup"]=abs(sup(u_h)-sup(u_ref))/abs(sup(u_ref))
109 m0_h=integrate(u_h_e)
110 m1_h=integrate(x*u_h_e)/m0_h
111 m2_h=integrate(length(x-m1_h)**2*u_h_e)
112
113 m0=integrate(u_ref)
114 m1=integrate(x*u_ref)/m0_h
115 m2=integrate(length(x-m1)**2*u_ref)
116
117 out["m0"]=abs(m0_h-m0)/abs(m0)
118 out["m1"]=length(m1_h-m1)/length(m1)
119 out["m2"]=abs(m2_h-m2)/abs(m2)
120
121 return out
122
123 #================
124
125 def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x", d_dir="x", a=1., CN=True):
126 """
127 dim - sparial dimension
128 s - width of initial profile
129 h - mesh size
130 """
131 v_c=c/a*getDirection(dim,c_dir)
132 v_d=d/a*getDirection(dim,d_dir)
133 v = (v_c+v_d)
134 E=b/a
135 if VERBOSITY:
136 print("="*100)
137 print("XX 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))
138 print("="*100)
139 print("initial width s = ",s)
140 print("diffusion = ",E)
141 print("total velocity = ",v)
142 print("tend = ", tend)
143 print("tolerance = ",TOL)
144 print("number of elements over s =",s/h)
145 b0=sqrt(- log(TAU) * 4*(s**2+E*tend))
146 b1=sqrt(- log(TAU)) * 2*s
147 X0_0=max(b1,-v[0]*tend + b0)
148 X0_1=max(b1,-v[1]*tend + b0)
149 l_0=X0_0+max(v[0]*tend + b0 , b1)
150 l_1=X0_1+max(v[1]*tend + b0 , b1)
151 NE_0=max(int(l_0/h+0.5),1)
152 NE_1=max(int(l_1/h+0.5),1)
153 if dim == 2:
154 if VERBOSITY: print("%d x %d grid over %e x %e with element size %e."%(NE_0,NE_1,l_0,l_1,h))
155 if NE_0*NE_1 > NE_MAX:
156 raise ValueError("too many elements %s."%(NE_0*NE_1,))
157 dom=Rectangle(n0=NE_0,n1=NE_1,l0=l_0,l1=l_1)
158 x0=[X0_0, X0_1]
159 else:
160 X0_2=max(b1,-v[2]*tend + b0)
161 l_2=X0_2+max(v[2]*tend + b0 , b1)
162 NE_2=max(int(l_2/h+0.5),1)
163 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))
164 if NE_0*NE_1*NE_2 > NE_MAX:
165 raise ValueError("too many elements %s."%(NE_0*NE_1*NE_2,))
166 dom=Brick(n0=NE_0,n1=NE_1, ne2=NE_2, l0=l_0,l1=l_1, l2=l_2)
167 x0=[X0_0, X0_1, X0_2]
168 if VERBOSITY:
169 print("initial location = ",x0)
170 print("XX", interpolate(uRef(dom,0.,E,s,v,x0), FunctionOnBoundary(dom)))
171
172 fc_BE=TransportPDE(dom,numEquations=1,useBackwardEuler=True)
173 fc_BE.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
174 fc_BE.getSolverOptions().setVerbosity(VERBOSITY)
175 fc_BE.getSolverOptions().setTolerance(TOL)
176 #
177 fc_BE.getSolverOptions().setPreconditioner(fc_BE.getSolverOptions().GAUSS_SEIDEL)
178 fc_BE.getSolverOptions().setNumSweeps(5)
179 if VERBOSITY: print("Backward Euler Transport created")
180
181 fc_CN=TransportPDE(dom,numEquations=1,useBackwardEuler=False)
182 fc_CN.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
183 fc_CN.getSolverOptions().setVerbosity(VERBOSITY)
184 fc_CN.getSolverOptions().setTolerance(TOL)
185
186 #fc_CN.getSolverOptions().setPreconditioner(fc_CN.getSolverOptions().GAUSS_SEIDEL)
187 fc_CN.getSolverOptions().setNumSweeps(2)
188 if VERBOSITY: print("Crank Nicolson Transport created")
189 dt_CN=fc_CN.getSafeTimeStepSize()
190 if VERBOSITY: print("time step size by Crank Nicolson=",dt_CN)
191
192 U0=uRef(dom,0,E,s,v,x0)
193 U0_e=uRef(dom,0,E,s,v,x0,True)
194 fc_CN.setInitialSolution(U0)
195 fc_BE.setInitialSolution(U0)
196 initial_error_L2=sqrt(integrate((U0-U0_e)**2))
197 if VERBOSITY:
198 print("initial Lsup = ",Lsup(U0), Lsup(U0_e))
199 print("initial integral = ",integrate(U0_e))
200 print("initial error = ",initial_error_L2)
201 print("used time step size =",dt)
202
203 if not CN:
204 n=int(ceil(tend/dt))
205 if VERBOSITY:
206 print("Solve Backward Euler:")
207 print("substeps : ",n)
208 t0=clock()
209 for i in range(n): u=fc_BE.getSolution(dt)
210 t0=clock()-t0
211 else:
212 if VERBOSITY: print("Solve Crank Nicolson:")
213 dt=dt_CN
214 t0=clock()
215 u=fc_CN.getSolution(tend)
216 t0=clock()-t0
217 out=QUALITY(u,uRef(dom,tend,E,s,v,x0,True))
218 print("XX", interpolate(uRef(dom,tend,E,s,v,x0), FunctionOnBoundary(dom)))
219 out['time']=t0
220 out['tend']=tend
221 out['dt']=dt
222 out['dx']=h
223 if abs(b)>0:
224 out["peclet"]=length(v)*s/b
225 else:
226 out["peclet"]=9999999.
227 # saveVTK("bb.vtu",u0=U0,u_CN=u_CN, uRef=uRef(dom,dt2,E,s,v,X0) )
228 return out
229
230 # (s, peclet, b, h0) -> error < 0.01
231 test_set = ( (0.05, 1., 1., 0.024), )
232 test_set = ( (0.05, 100000., 1., 0.024), )
233
234 if False:
235 S_MAX=0.5/sqrt(-log(TAU))/2
236 s=0.05
237 peclet = 1000.
238 b=1.
239 c=peclet*b/s
240 h=0.1/4*1.2/1.25
241 dt=6.250000e-10
242
243 print(XXX(DIM,dt,dt,s=s,h=h,b=a*b,c=a*c,d=0,c_dir="x", d_dir="x", CN=True))
244 1/0
245
246 for tst in test_set:
247 s=tst[0]
248 peclet=tst[1]
249 b=tst[2]
250 h0=tst[3]
251 c=peclet*b/s
252
253 # find appropraiate tend:
254 result=XXX(DIM,1e-99,1.,s=s,h=h0,b=a*b,c=a*c,d=0,c_dir="x", d_dir="x", CN=True)
255 tend=result['dt']
256
257 f_test=[ 1 , 2, 4 ]
258 f_test=[ 1, 2, 4, 8 ]
259 out=""
260 tab_name="dt"
261 tab_name="tend"
262 tab_name="error"
263 dt_s=[]
264 for f_h in f_test:
265 out+="h0/%s "%f_h
266 h=h0/f_h
267 result=XXX(DIM,tend,tend,s=s,h=h,b=a*b,c=a*c,d=0,c_dir="x", d_dir="x", CN=True)
268 out+=", %e"%result[tab_name]
269 print("XX",result)
270 dt_s.insert(0,result['dt'])
271 for i in range(len(f_test)-len(dt_s)): out+=", "
272 for dt in dt_s:
273 result=XXX(DIM,tend,dt,s=s,h=h,b=a*b,c=a*c,d=0,c_dir="x", d_dir="x", CN=False)
274 print("XX",result)
275 out+=", %e"%result[tab_name]
276 out+="\n"
277 header="h\dt , "
278 for dt in dt_s: header+=", %e"%dt
279 out=header+"\n"+out
280 print(out)

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26