# Diff of /trunk/finley/test/python/FCT_benchmark.py

revision 3070 by gross, Wed Jun 2 10:17:00 2010 UTC revision 3071 by gross, Wed Jul 21 05:37:30 2010 UTC
# Line 25  __url__="https://launchpad.net/escript-f Line 25  __url__="https://launchpad.net/escript-f
25  #  #
26  #     we solve a* U_{,t} - b *u_{,ii} + c_i u_{,i} + (d_i* u)_{,i}}=0  #     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) )  #               U(0)= U0 * exp ( - |x-x_0(t)|^2/(4*s**2) )
29  #  #
30  #  with a>0, b>=0, s>0  #  with a>0, b>=0, s>0
31  #  #
32  #  we set E=b/a v=c/a w=d/a  #  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)) )  #  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  #  #
36  #                 with x_0(t) = X0 + (v+w)*t  #                 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  #    the region |x-x_0(t)|^2/(4*(s**2+E*t)) < - log(TAU)  is within the domain for all time
# domain where TAU is a given relative tolerance
40  #  #
41  #  #
42  #    this holds if  #    this holds if
43  #  #
44  #       |x_i-X0_i-(v_i+w_i)*tend | < - log(TAU) * 4*(s+E*tend)  #       |x_i-X0_i-(v_i+w_i)*tend | < sqrt(- log(TAU) * 4*(s**2+E*tend))=b0  and
45  #       |x_i-X0_i | < - log(TAU) * 4*s  #       |x_i-X0_i | < sqrt(- log(TAU)) * 2*s = b1 implies 0<=x_i<=l_i
46  #  #
47  #        X0_i=min(- log(TAU) * 4*s,-(v_i+w_i)*tend - log(TAU) * 4*(s+E*tend))  from math import pi, ceil
48  #        l_i=X0_i+min((v_i+w_i)*tend - log(TAU) * 4*(s+E*tend)) , - log(TAU) * 4*s)  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  #    #
53  DIM=2  DIM=2
54  NE_MAX=35000  NE_MAX=300000
55  VERBOSITY=True  VERBOSITY=True
56  TOL=1.e-8  TOL=1.e-8
57  TAU=1e-10  TAU=1e-10
58  VTK_DIR="output"  VTK_DIR="output"
59  S_RAISE=1.5
60    #==================
61  S_MIN=0  S_MIN=0
S_MAX=0
62  TABS = [ 'dx', 'dt', 'peclet', 'error',  'sup', 'integral',  'center', 'width', 'time' ]  TABS = [ 'dx', 'dt', 'peclet', 'error',  'sup', 'integral',  'center', 'width', 'time' ]
63  from math import pi, ceil
64  from time import time as clock  a=1.
65  from esys.finley import Rectangle, Brick  #==================
from esys.escript import *
from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
S_MAX=-0.5/log(TAU)/4
66
67
68  mkDir(VTK_DIR)  mkDir(VTK_DIR)
# Line 74  def uRef(dom,t,E,s,v,x0, onElements=Fals Line 73  def uRef(dom,t,E,s,v,x0, onElements=Fals
73      else:      else:
74        x=dom.getX()        x=dom.getX()
75      X=x0[:dom.getDim()]+v[:dom.getDim()]*t      X=x0[:dom.getDim()]+v[:dom.getDim()]*t
76      u=(s/(s+E*t))**(dom.getDim()/2.) * exp(-length(x-X)**2/(4*(s+E*t)))      u=(s**2/(s**2+E*t))**(dom.getDim()/2.) * exp(-length(x-X)**2/(4*(s**2+E*t)))
77      return u      return u
78
NE_test = [ 10, 20, 40, 80, 160 ]
NE_test = [ 10, 20, 40, 80, 160, 320]
NE_test = [ 320 *2 ]
b_test=[ 0., 0.001, 0.1, 1., 10, 100.]
b_test=[ 1.0 *0]
c_test=[ 0., 0.001, 0.1, 1., 10, 100. ]
c_test=[ 100.]
d_test=[ 0.0  ]
s_test=[0.0005 ]
# s=[ 0.001, 0.003, 0.01 ]

a=1.
X0=[0.5,0.5,0.5]
step=0.1
s_max=0.01  # don't use an s which is bigger than this as the profile will not fit into the domian

79
80  def getDirection(dim, d="x"):  def getDirection(dim, d="x"):
81       k=kronecker(dim)       k=kronecker(dim)
# Line 113  def getDirection(dim, d="x"): Line 96  def getDirection(dim, d="x"):
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       dom=u_h.getDomain()       u_h_e=interpolate(u_h,u_ref.getFunctionSpace())
101       x=dom.getX()       x=u_ref.getFunctionSpace().getX()
102       out = {}       out = {}
103         out["error"]=sqrt(integrate((u_ref-u_h_e)**2))/sqrt(integrate(u_ref**2))
104       out["sup"]=abs(sup(u_h)-sup(u_ref))/abs(sup(u_ref))       out["sup"]=abs(sup(u_h)-sup(u_ref))/abs(sup(u_ref))
105       i_h=integrate(u_h)       m0_h=integrate(u_h_e)
106       i_ref=integrate(u_ref)       m1_h=integrate(x*u_h_e)/m0_h
107       out["integral"]=abs(i_h-i_ref)/abs(i_ref)       m2_h=integrate(length(x-m1_h)**2*u_h_e)
108       x0h=integrate(x*u_h)/i_h
109       x0ref=integrate(x*u_ref)/i_ref       m0=integrate(u_ref)
110       out["center"]=length(x0h-x0ref)/(1.+step)       m1=integrate(x*u_ref)/m0_h
111         m2=integrate(length(x-m1)**2*u_ref)
112       sigma_h=sqrt(integrate(length(x-x0h)**2*u_h))
113       sigma_ref=sqrt(integrate(length(x-x0ref)**2*u_ref))       out["m0"]=abs(m0_h-m0)/abs(m0)
114       out["width"]=abs(sigma_h-sigma_ref)/abs(sigma_ref)       out["m1"]=length(m1_h-m1)/length(m1)
115       out["error"]=sqrt(integrate((u_ref-u_h)**2))       out["m2"]=abs(m2_h-m2)/abs(m2)
116
117       return out       return out
118
119    #================
120
121    def XXX(dim,tend,dt, s, h,b,c,d,c_dir="x", d_dir="x", a=1., CN=True):
122  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir="x", a=1.):      """
123         dim - sparial dimension
124         s - width of initial profile
125         h - mesh size
126        """
127      v_c=c/a*getDirection(dim,c_dir)      v_c=c/a*getDirection(dim,c_dir)
128      v_d=d/a*getDirection(dim,d_dir)      v_d=d/a*getDirection(dim,d_dir)
129      v = (v_c+v_d)      v = (v_c+v_d)
130      E=b/a      E=b/a
if abs(E)>0:
if length(v)>0:
tend=min(s/E, s/length(v))
else:
tend=s/E
else:
if length(v)>0:
tend=s/length(v)
else:
tend=0
tend*=(1.-S_RAISE)
131      if VERBOSITY:      if VERBOSITY:
132             print "="*100             print "="*100
133             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)             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 "XX number of elements over s =",s/h             print "number of elements over s =",s/h
141      X0_0=min(- log(TAU) * 4*s,-v[0]*tend - log(TAU) * 4*(s+E*tend))      b0=sqrt(- log(TAU) * 4*(s**2+E*tend))
142      X0_1=min(- log(TAU) * 4*s,-v[1]*tend - log(TAU) * 4*(s+E*tend))      b1=sqrt(- log(TAU)) * 2*s
143      l_0=X0_0+min(v[0]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)      X0_0=max(b1,-v[0]*tend + b0)
144      l_1=X0_1+min(v[1]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)      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      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:
# Line 174  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir Line 153  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir
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=min(- log(TAU) * 4*s,-v[2]*tend - log(TAU) * 4*(s+E*tend))         X0_2=max(b1,-v[2]*tend + b0)
157         l_2=X0_2+min(v[2]*tend - log(TAU) * 4*(s+E*tend) , - log(TAU) * 4*s)         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:
# Line 184  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir Line 163  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir
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))
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)
170      fc_BE.getSolverOptions().setVerbosity(VERBOSITY)      fc_BE.getSolverOptions().setVerbosity(VERBOSITY)
171      fc_BE.getSolverOptions().setTolerance(TOL)      fc_BE.getSolverOptions().setTolerance(TOL)
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(1)        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)
# Line 204  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir Line 184  def XXX(dim, s, h,b,c,d,c_dir="x", d_dir
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      if tend>0:
188         dt=main(dt_CN, tend)      U0=uRef(dom,0,E,s,v,x0)
189      else:      U0_e=uRef(dom,0,E,s,v,x0,True)
dt=1.
U0=uRef(dom,0,E,s,v,X0)
U0_e=uRef(dom,0,E,s,v,X0,True)
190      fc_CN.setInitialSolution(U0)      fc_CN.setInitialSolution(U0)
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 "used time step size =",dt        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
198
199            if not CN:
200           n=int(ceil(tend/dt))
201      return         if VERBOSITY:
202      if False:            print "Solve Backward Euler:"
203      # dt2=min(dt2,20*dt)            print "substeps : ",n
204      if VERBOSITY: print "Solve Backward Euler:"         t0=clock()
205          t0=clock()         for i in range(n): u=fc_BE.getSolution(dt)
206          u_BE=fc_BE.getSolution(10.*dt)         t0=clock()-t0
207          u_BE=fc_BE.getSolution(10.*dt)      else:
208      t0=clock()-t0         if VERBOSITY: print "Solve Crank Nicolson:"
209      out=QUALITY(u_BE,uRef(dom,dt2,E,s,v,X0))         dt=dt_CN
210      #             t0=clock()
211      # print integrate(U0), U0             u=fc_CN.getSolution(tend)
212      # print integrate(u_BE), u_BE         t0=clock()-t0
213      # print integrate(uRef(dom,dt2,E,s,v,X0)), uRef(dom,dt2,E,s,v,X0)      out=QUALITY(u,uRef(dom,tend,E,s,v,x0,True))
214      # saveVTK("bb_be.vtu",u0=U0,u_BE=u_BE, uRef=uRef(dom,dt2,E,s,v,X0) )      print "XX", interpolate(uRef(dom,tend,E,s,v,x0), FunctionOnBoundary(dom))
215      # 1/0      out['time']=t0
216      #      out['tend']=tend
else:
if VERBOSITY: print "Solve Crank Nicolson:"
t0=clock()
u_CN=fc_CN.getSolution(dt)
t0=clock()-t0
out=QUALITY(u_CN,uRef(dom,dt,E,s,v,X0))
out['time']=t0/NN
217      out['dt']=dt      out['dt']=dt
218      out['dx']=1./NE      out['dx']=h
219      if abs(b)>0:      if abs(b)>0:
220         out["peclet"]=length(v)/(b*NE)         out["peclet"]=length(v)*s/b
221      else:      else:
222          out["peclet"]=9999999.          out["peclet"]=9999999.
223      # saveVTK("bb.vtu",u0=U0,u_CN=u_CN, uRef=uRef(dom,dt2,E,s,v,X0) )      # saveVTK("bb.vtu",u0=U0,u_CN=u_CN, uRef=uRef(dom,dt2,E,s,v,X0) )
224      return out      return out
225
226    # (s, peclet, b, h0)  -> error < 0.01
227  for s_fac in [ 1., 0.5, 0.25, 0.125 ] :  test_set = ( (0.05, 1., 1., 0.024),  )
228     for h_fac in [1., 0.5, 0.25, 0.125 ]:  test_set = ( (0.05, 100000., 1., 0.024),  )
229        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")
230  1/0  if False:
231        S_MAX=0.5/sqrt(-log(TAU))/2
232  txt=""      s=0.05
233  for t in TABS: txt+=" "+t      peclet = 1000.
234  for b in b_test:      b=1.
235    for c in c_test:      c=peclet*b/s
236      for d in d_test:      h=0.1/4*1.2/1.25
237        for s in s_test:      dt=6.250000e-10
238      for NE in NE_test:
239          dat=XXX(NE, a,b,c,d,c_dir="x", d_dir="x",s=s)      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              txt+="\n"      1/0
241              for t in TABS: txt+=" %3.2e"%dat[t]
242                for tst in test_set:
243  print txt       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