/[escript]/trunk/finley/test/python/FCT_test2.py
ViewVC logotype

Annotation of /trunk/finley/test/python/FCT_test2.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2344 - (hide annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 4166 byte(s)
Change __url__ to launchpad site

1 ksteube 1811
2     ########################################################
3 gross 1400 #
4 ksteube 1811 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-2008 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     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1811
22     #
23 gross 1400 # upwinding test moving a Gaussian hill around
24     #
25     # we solve U_,t - E *u_,ii + v_i u_,i =0 (E is small)
26     #
27     # the solution is given as u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) )
28     #
29     # where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and
30     #
31     # x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3
32     #
33     # the solution is started from some time T0>0.
34     #
35     # We are using five quality messurements for u_h
36     #
37     # - inf(u_h) > 0
38     # - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1
39     # - integrate(u_h) ~ 1
40     # - | x_0h-x_0 | ~ 0 where x_0h = integrate(x*u_h)
41     # - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 )
42     #
43     #
44 gross 1410 from esys.escript import *
45     from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
46     from esys.finley import Rectangle, Brick
47 gross 1400 from math import pi, ceil
48 gross 1410 NE=128
49 gross 1476 NE=64
50 gross 1400 DIM=2
51     THETA=0.5
52     OMEGA0=1.
53     ALPHA=pi/4
54     T0=0.5*pi
55 gross 1410 T_END=2.5*pi
56     dt=1e-3*10
57     E=1.e-3
58 gross 1476 TEST_SUPG=False or True
59 gross 1400
60 gross 1410
61 gross 1400 def getCenter(t):
62     if DIM==2:
63     x0=[cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5]
64 gross 1410 x0=[-sin(OMEGA0*t)*0.5,cos(OMEGA0*t)*0.5]
65 gross 1400 else:
66     x0=[cos(ALPHA)*cos(OMEGA0*t)*0.5,-sin(OMEGA0*t)*0.5,-sin(ALPHA)*cos(OMEGA0*t)*0.5]
67     return x0
68     def QUALITY(t,u_h):
69     dom=u_h.getDomain()
70     x=dom.getX()
71     a=inf(u_h)
72     b=sup(u_h)*(4*pi*E*t)**(DIM/2.)-1.
73     c=integrate(u_h,Function(dom))-1.
74     x0=getCenter(t)
75     x0h=integrate(x*u_h,Function(dom))
76     d=length(x0-x0h)
77     sigma_h2=sqrt(integrate(length(x-x0h)**2 * u_h, Function(dom)))
78     if DIM == 3: sigma_h2*=sqrt(2./3.)
79     e=sigma_h2/sqrt(4*E*t)-1
80 gross 1410 # return a,b,c,e,1./(4*pi*E*t)**(DIM/2.)
81 gross 1476 return d,e
82 gross 1410 # return a,b,c,d,e
83 gross 1400
84    
85 gross 1368
86    
87 gross 1400 if DIM==2:
88     dom=Rectangle(NE,NE)
89     else:
90     dom=Brick(NE,NE,NE)
91     dom.setX(2*dom.getX()-1)
92    
93     # set initial value
94 gross 1370 x=dom.getX()
95 gross 1400 u0=1/(4.*pi*E*T0)**(DIM/2.)*exp(-length(dom.getX()-getCenter(T0))**2/(4.*E*T0))
96 gross 1371
97 gross 1400 print "QUALITY ",QUALITY(T0,u0)
98    
99     x=Function(dom).getX()
100     if DIM == 2:
101 gross 1410 V=OMEGA0*(x[0]*[0,-1]+x[1]*[1,0])
102 gross 1400 else:
103     V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
104     #===================
105     fc=TransportPDE(dom,num_equations=1,theta=THETA)
106     x=Function(dom).getX()
107     fc.setValue(M=Scalar(1.,Function(dom)),C=V,A=-Scalar(E,Function(dom))*kronecker(dom))
108     #==============
109     if TEST_SUPG:
110     supg=LinearSinglePDE(dom)
111     supg.setValue(D=1.)
112     supg.setSolverMethod(supg.LUMPING)
113     dt_supg=1./(1./inf(dom.getSize()/length(V))+1./inf(dom.getSize()**2/E))*0.3
114     u_supg=u0*1.
115    
116 gross 1370 c=0
117 gross 1400 saveVTK("u.%s.xml"%c,u=u0)
118     fc.setInitialSolution(u0)
119     t=T0
120 gross 1410 while t<T_END:
121 gross 1370 print "time step t=",t+dt
122     u=fc.solve(dt)
123 gross 1400 if TEST_SUPG:
124     #========== supg tests ================
125     nn=max(ceil(dt/dt_supg),1.)
126     dt2=dt/nn
127     nnn=0
128     while nnn<nn :
129     supg.setValue(X=-dt2/2*E*grad(u_supg),Y=u_supg+dt2/2*inner(V,grad(u_supg)))
130     u2=supg.getSolution()
131     supg.setValue(X=-dt2*E*grad(u2),Y=u_supg+dt2*inner(V,grad(u2)))
132     u_supg=supg.getSolution()
133     nnn+=1
134 gross 1370 c+=1
135     t+=dt
136 gross 1476 print "QUALITY FCT: time = %s pi"%(t/pi),QUALITY(t,u),
137 gross 1400 if TEST_SUPG:
138     print "QUALITY SUPG: ",QUALITY(t,u_supg)
139     # saveVTK("u.%s.xml"%c,u=u,u_supg=u_supg)
140     else:
141     # saveVTK("u.%s.xml"%c,u=u)
142     pass
143     # if c == 20: 1/0

  ViewVC Help
Powered by ViewVC 1.1.26