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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (22 months ago) by jfenwick
File MIME type: text/x-python
File size: 8789 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26