# Annotation of /branches/symbolic_from_3470/dudley/test/python/FCT_benchmark.py

Revision 3789 - (hide annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 8620 byte(s)
```Fast forward to latest trunk revision 3788.

```
 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 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 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 jfenwick 3087 from esys.dudley import Rectangle, Brick 50 gross 3071 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 caltinay 3789 raise ValueError("Cannot identify direction %s"%d) 98 gross 3033 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 caltinay 3789 print("="*100) 133 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 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 print("number of elements over s =",s/h) 141 gross 3071 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 caltinay 3789 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 gross 3033 if NE_0*NE_1 > NE_MAX: 152 caltinay 3789 raise ValueError("too many elements %s."%(NE_0*NE_1,)) 153 gross 3033 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 caltinay 3789 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 gross 3033 if NE_0*NE_1*NE_2 > NE_MAX: 161 caltinay 3789 raise ValueError("too many elements %s."%(NE_0*NE_1*NE_2,)) 162 gross 3033 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 caltinay 3789 print("initial location = ",x0) 166 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 caltinay 3789 if VERBOSITY: print("Backward Euler Transport created") 176 gross 3033 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 caltinay 3789 if VERBOSITY: print("Crank Nicolson Transport created") 185 gross 3033 dt_CN=fc_CN.getSafeTimeStepSize() 186 caltinay 3789 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 caltinay 3789 print("initial Lsup = ",Lsup(U0), Lsup(U0_e)) 195 print("initial integral = ",integrate(U0_e)) 196 print("initial error = ",initial_error_L2) 197 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 caltinay 3789 print("Solve Backward Euler:") 203 print("substeps : ",n) 204 gross 3071 t0=clock() 205 for i in range(n): u=fc_BE.getSolution(dt) 206 t0=clock()-t0 207 gross 3033 else: 208 caltinay 3789 if VERBOSITY: print("Solve Crank Nicolson:") 209 gross 3071 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 caltinay 3789 print("XX", interpolate(uRef(dom,tend,E,s,v,x0), FunctionOnBoundary(dom))) 215 gross 3071 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 caltinay 3789 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 gross 3071 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 caltinay 3789 print("XX",result) 266 gross 3071 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 caltinay 3789 print("XX",result) 271 gross 3071 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 caltinay 3789 print(out)

Name Value
svn:executable *