Diff of /branches/symbolic_from_3470/dudley/test/python/FCT_benchmark.py

revision 3788 by caltinay, Tue Mar 15 04:23:54 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 94  def getDirection(dim, d="x"): Line 94  def getDirection(dim, d="x"):
94       elif d=="xyz" and dim>2:       elif d=="xyz" and dim>2:
95           return (k[1]+k[2]+k[0])/sqrt(3.)           return (k[1]+k[2]+k[0])/sqrt(3.)
96       else:       else:
97           raise ValueError,"Cannot identify direction %s"%d           raise ValueError("Cannot identify direction %s"%d)
98
99  def QUALITY(u_h,u_ref):  def QUALITY(u_h,u_ref):
100       u_h_e=interpolate(u_h,u_ref.getFunctionSpace())       u_h_e=interpolate(u_h,u_ref.getFunctionSpace())
# Line 129  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x Line 129  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x
129      v = (v_c+v_d)      v = (v_c+v_d)
130      E=b/a      E=b/a
131      if VERBOSITY:      if VERBOSITY:
132             print "="*100             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)             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             print("="*100)
135             print "initial width s = ",s             print("initial width s = ",s)
136             print "diffusion = ",E             print("diffusion = ",E)
137             print "total velocity = ",v             print("total velocity = ",v)
138             print "tend = ", tend             print("tend = ", tend)
139             print "tolerance = ",TOL             print("tolerance = ",TOL)
140             print "number of elements over s =",s/h             print("number of elements over s =",s/h)
141      b0=sqrt(- log(TAU) * 4*(s**2+E*tend))      b0=sqrt(- log(TAU) * 4*(s**2+E*tend))
142      b1=sqrt(- log(TAU)) * 2*s      b1=sqrt(- log(TAU)) * 2*s
143      X0_0=max(b1,-v[0]*tend + b0)      X0_0=max(b1,-v[0]*tend + b0)
# Line 147  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x Line 147  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x
147      NE_0=max(int(l_0/h+0.5),1)      NE_0=max(int(l_0/h+0.5),1)
148      NE_1=max(int(l_1/h+0.5),1)      NE_1=max(int(l_1/h+0.5),1)
149      if dim == 2:      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)          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:          if NE_0*NE_1 > NE_MAX:
152         raise ValueError,"too many elements %s."%(NE_0*NE_1,)         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)          dom=Rectangle(n0=NE_0,n1=NE_1,l0=l_0,l1=l_1)
154          x0=[X0_0, X0_1]          x0=[X0_0, X0_1]
155      else:      else:
156         X0_2=max(b1,-v[2]*tend + b0)         X0_2=max(b1,-v[2]*tend + b0)
157         l_2=X0_2+max(v[2]*tend + b0 , b1)         l_2=X0_2+max(v[2]*tend + b0 , b1)
158         NE_2=max(int(l_2/h+0.5),1)         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)         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:         if NE_0*NE_1*NE_2 > NE_MAX:
161        raise ValueError,"too many elements %s."%(NE_0*NE_1*NE_2,)        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)         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]         x0=[X0_0, X0_1, X0_2]
164      if VERBOSITY:      if VERBOSITY:
165          print "initial location = ",x0          print("initial location = ",x0)
166      print "XX", interpolate(uRef(dom,0.,E,s,v,x0), FunctionOnBoundary(dom))      print("XX", interpolate(uRef(dom,0.,E,s,v,x0), FunctionOnBoundary(dom)))
167
168      fc_BE=TransportPDE(dom,numEquations=1,useBackwardEuler=True)      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)      fc_BE.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
# Line 172  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x Line 172  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x
172      #      #
173      fc_BE.getSolverOptions().setPreconditioner(fc_BE.getSolverOptions().GAUSS_SEIDEL)      fc_BE.getSolverOptions().setPreconditioner(fc_BE.getSolverOptions().GAUSS_SEIDEL)
174      fc_BE.getSolverOptions().setNumSweeps(5)        fc_BE.getSolverOptions().setNumSweeps(5)
175      if VERBOSITY: print "Backward Euler Transport created"      if VERBOSITY: print("Backward Euler Transport created")
176
177      fc_CN=TransportPDE(dom,numEquations=1,useBackwardEuler=False)      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)      fc_CN.setValue(M=a, A=-b*kronecker(dom), B=-v_d*a, C=-v_c*a)
# Line 181  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x Line 181  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x
181
182      #fc_CN.getSolverOptions().setPreconditioner(fc_CN.getSolverOptions().GAUSS_SEIDEL)      #fc_CN.getSolverOptions().setPreconditioner(fc_CN.getSolverOptions().GAUSS_SEIDEL)
183      fc_CN.getSolverOptions().setNumSweeps(2)        fc_CN.getSolverOptions().setNumSweeps(2)
184      if VERBOSITY: print "Crank Nicolson Transport created"      if VERBOSITY: print("Crank Nicolson Transport created")
185      dt_CN=fc_CN.getSafeTimeStepSize()      dt_CN=fc_CN.getSafeTimeStepSize()
186      if VERBOSITY: print "time step size by Crank Nicolson=",dt_CN      if VERBOSITY: print("time step size by Crank Nicolson=",dt_CN)
187
188      U0=uRef(dom,0,E,s,v,x0)      U0=uRef(dom,0,E,s,v,x0)
189      U0_e=uRef(dom,0,E,s,v,x0,True)      U0_e=uRef(dom,0,E,s,v,x0,True)
# Line 191  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x Line 191  def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x
191      fc_BE.setInitialSolution(U0)      fc_BE.setInitialSolution(U0)
192      initial_error_L2=sqrt(integrate((U0-U0_e)**2))      initial_error_L2=sqrt(integrate((U0-U0_e)**2))
193      if VERBOSITY:      if VERBOSITY:
194        print "initial Lsup = ",Lsup(U0), Lsup(U0_e)        print("initial Lsup = ",Lsup(U0), Lsup(U0_e))
195        print "initial integral = ",integrate(U0_e)        print("initial integral = ",integrate(U0_e))
196        print "initial error = ",initial_error_L2        print("initial error = ",initial_error_L2)
197        print "used time step size =",dt        print("used time step size =",dt)
198
199      if not CN:      if not CN:
200         n=int(ceil(tend/dt))         n=int(ceil(tend/dt))
201         if VERBOSITY:         if VERBOSITY:
202            print "Solve Backward Euler:"            print("Solve Backward Euler:")
203            print "substeps : ",n            print("substeps : ",n)
204         t0=clock()         t0=clock()
205         for i in range(n): u=fc_BE.getSolution(dt)         for i in range(n): u=fc_BE.getSolution(dt)
206         t0=clock()-t0         t0=clock()-t0
207      else:      else:
208         if VERBOSITY: print "Solve Crank Nicolson:"         if VERBOSITY: print("Solve Crank Nicolson:")
209         dt=dt_CN         dt=dt_CN
210         t0=clock()         t0=clock()
211         u=fc_CN.getSolution(tend)         u=fc_CN.getSolution(tend)
212         t0=clock()-t0         t0=clock()-t0
213      out=QUALITY(u,uRef(dom,tend,E,s,v,x0,True))      out=QUALITY(u,uRef(dom,tend,E,s,v,x0,True))
214      print "XX", interpolate(uRef(dom,tend,E,s,v,x0), FunctionOnBoundary(dom))      print("XX", interpolate(uRef(dom,tend,E,s,v,x0), FunctionOnBoundary(dom)))
215      out['time']=t0      out['time']=t0
216      out['tend']=tend      out['tend']=tend
217      out['dt']=dt      out['dt']=dt
# Line 236  if False: Line 236  if False:
236      h=0.1/4*1.2/1.25      h=0.1/4*1.2/1.25
237      dt=6.250000e-10      dt=6.250000e-10
238
239      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)      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      1/0
241
242  for tst in test_set:  for tst in test_set:
# Line 262  for tst in test_set: Line 262  for tst in test_set:
262          h=h0/f_h          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)          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]          out+=", %e"%result[tab_name]
265          print "XX",result          print("XX",result)
266          dt_s.insert(0,result['dt'])          dt_s.insert(0,result['dt'])
267          for i in range(len(f_test)-len(dt_s)): out+=", "          for i in range(len(f_test)-len(dt_s)): out+=", "
268          for dt in dt_s:          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)              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              print("XX",result)
271              out+=", %e"%result[tab_name]              out+=", %e"%result[tab_name]
272          out+="\n"          out+="\n"