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

trunk/finley/test/python/tp.py revision 1393 by gross, Mon Jan 21 06:20:54 2008 UTC trunk/finley/test/python/FCT_test2.py revision 1400 by gross, Thu Jan 24 06:04:31 2008 UTC
# Line 1  Line 1
1    #
2    #  upwinding test moving a Gaussian hill around
3    #
4    #     we solve U_,t - E *u_,ii + v_i u_,i =0 (E is small)
5    #
6    #  the solution is given as   u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) )
7    #
8    #   where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and
9    #
10    #         x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3
11    #
12    #  the solution is started from some time T0>0.
13    #
14    #  We are using five quality messurements for u_h
15    #
16    #     - inf(u_h) > 0
17    #     - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1
18    #     - integrate(u_h) ~ 1
19    #     - | x_0h-x_0 | ~ 0    where x_0h = integrate(x*u_h)
20    #     - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )
21    #
22    #
23    from math import pi, ceil
24    NE=128/2
25    DIM=2
26    THETA=0.5
27    OMEGA0=1.
28    ALPHA=pi/4
29    T0=0.5*pi
30    E=5e-4 * 2
31    TEST_SUPG=False
32
33    def getCenter(t):
34       if DIM==2:
35          x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5]
36       else:
37          x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5]
38       return x0
39    def QUALITY(t,u_h):
40         dom=u_h.getDomain()
41         x=dom.getX()
42         a=inf(u_h)
43         b=sup(u_h)*(4*pi*E*t)**(DIM/2.)-1.
44         c=integrate(u_h,Function(dom))-1.
45         x0=getCenter(t)
46         x0h=integrate(x*u_h,Function(dom))
47         d=length(x0-x0h)
48         sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom)))
49         if DIM == 3: sigma_h2*=sqrt(2./3.)
50         e=sigma_h2/sqrt(4*E*t)-1
51         return a,b,c,d,e
52
53
54  from esys.escript import *  from esys.escript import *
55    from esys.escript.linearPDEs import LinearSinglePDE
56
57  class TransportPDE(object):  class TransportPDE(object):
58       def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):       def __init__(self,domain,num_equations=1,theta=0.,dt_max=-1.,trace=True):
# Line 56  class TransportPDE(object): Line 110  class TransportPDE(object):
110       def setInitialSolution(self,u):       def setInitialSolution(self,u):
111           self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))           self.__transport_problem.setInitialValue(interpolate(u,self.getFunctionSpace()))
112       def solve(self,dt):       def solve(self,dt):
113         return self.__transport_problem.solve(self.__source,dt,{"verbose" : True })         return self.__transport_problem.solve(self.__source,dt,{"verbose" : True , "tolerance" : 1.e-8})
114  from esys.finley import Rectangle
115
116    from esys.finley import Rectangle, Brick
117
118
119  dom=Rectangle(40,20,l0=2)  if DIM==2:
120  fc=TransportPDE(dom,num_equations=1,theta=1.0,dt_max=2.5e-2/10)    dom=Rectangle(NE,NE)
121  fc.setValue(M=Scalar(1.,Function(dom)),C=Scalar(1.,Function(dom))*[-1.,0])  else:
122      dom=Brick(NE,NE,NE)
123    dom.setX(2*dom.getX()-1)
124
125    # set initial value
126  x=dom.getX()  x=dom.getX()
127  x_0=[0.3,0.3]  u0=1/(4.*pi*E*T0)**(DIM/2.)*exp(-length(dom.getX()-getCenter(T0))**2/(4.*E*T0))
sigma=0.08
u=1.
for i in range(dom.getDim()):
u=u*exp(-(x[i]-x_0[i])**2/sigma**2)
u/=Lsup(u)
128
129  u=whereNonPositive(abs(x[0]-0.4)-0.2)*whereNonPositive(abs(x[1]-0.5)-0.2)  print "QUALITY ",QUALITY(T0,u0)
c=0
saveVTK("u.%s.xml"%c,u=u)
fc.setInitialSolution(u)
130
131  dt=2.5e-2  x=Function(dom).getX()
132  t=0.  if DIM == 2:
133  while t<25*dt:     V=OMEGA0*(x[0]*[0,1]+x[1]*[-1,0])
134    else:
135       V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
136    #===================
137    fc=TransportPDE(dom,num_equations=1,theta=THETA)
138    x=Function(dom).getX()
139    fc.setValue(M=Scalar(1.,Function(dom)),C=V,A=-Scalar(E,Function(dom))*kronecker(dom))
140    #==============
141    if TEST_SUPG:
142       supg=LinearSinglePDE(dom)
143       supg.setValue(D=1.)
144       supg.setSolverMethod(supg.LUMPING)
145       dt_supg=1./(1./inf(dom.getSize()/length(V))+1./inf(dom.getSize()**2/E))*0.3
146       u_supg=u0*1.
147
148    c=0
149    saveVTK("u.%s.xml"%c,u=u0)
150    fc.setInitialSolution(u0)
151    dt=1.266539e-02*10
152    t=T0
153    while t<T0+10:
154      print "time step t=",t+dt        print "time step t=",t+dt
155      u=fc.solve(dt)        u=fc.solve(dt)
156      print "range u",inf(u),sup(u),integrate(u,Function(dom))      if TEST_SUPG:
157            #========== supg tests ================
158            nn=max(ceil(dt/dt_supg),1.)
159            dt2=dt/nn
160            nnn=0
161            while nnn<nn :
162                supg.setValue(X=-dt2/2*E*grad(u_supg),Y=u_supg+dt2/2*inner(V,grad(u_supg)))
163                u2=supg.getSolution()
164                supg.setValue(X=-dt2*E*grad(u2),Y=u_supg+dt2*inner(V,grad(u2)))
165                u_supg=supg.getSolution()
166                nnn+=1
167      c+=1      c+=1
saveVTK("u.%s.xml"%c,u=u)
168      t+=dt      t+=dt
169        print "QUALITY FCT: ",QUALITY(t,u)
170        if TEST_SUPG:
171           print "QUALITY SUPG: ",QUALITY(t,u_supg)
172           # saveVTK("u.%s.xml"%c,u=u,u_supg=u_supg)
173        else:
174           # saveVTK("u.%s.xml"%c,u=u)
175           pass
176        # if c == 20: 1/0

Legend:
 Removed from v.1393 changed lines Added in v.1400

 ViewVC Help Powered by ViewVC 1.1.26