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

Revision 3033 - (show annotations)
Wed Jun 2 10:17:00 2010 UTC (11 years 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 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 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 20 __url__= 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

Name Value
svn:executable *