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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26