← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/FCT_test1.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/FCT_test1.py
 
1
 
 
2
########################################################
 
3
#
 
4
# Copyright (c) 2003-2010 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-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
http://www.opensource.org/licenses/osl-3.0.php"""
 
20
__url__="https://launchpad.net/escript-finley"
 
21
 
 
22
#
 
23
#  upwinding test moving a Gaussian hill around 
 
24
#
 
25
#     we solve U_,t + v_i u_,i =0 
 
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
from esys.escript import *
 
45
from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE
 
46
from esys.dudley import Rectangle, Brick
 
47
from math import pi, ceil
 
48
NE=128
 
49
NE=4
 
50
DIM=2
 
51
THETA=0.5
 
52
OMEGA0=1.
 
53
ALPHA=pi/4
 
54
T0=0
 
55
T_END=2.*pi
 
56
dt=1e-3*10*10
 
57
E=1.e-3
 
58
TEST_SUPG=False or True
 
59
 
 
60
 
 
61
if DIM==2:
 
62
  dom=Rectangle(NE,NE)
 
63
else:
 
64
  dom=Brick(NE,NE,NE)
 
65
u0=dom.getX()[0]
 
66
# saveVTK("u.%s.vtu"%0,u=u0)
 
67
# print "XX"*80
 
68
dom.setX(2*dom.getX()-1)
 
69
 
 
70
# set initial value 
 
71
x=dom.getX()
 
72
r=sqrt(x[0]**2+(x[1]-1./3.)**2)
 
73
# u0=whereNegative(r-1./3.)*wherePositive(wherePositive(abs(x[0])-0.05)+wherePositive(x[1]-0.5))
 
74
 
 
75
x=Function(dom).getX()
 
76
if DIM == 2:
 
77
   V=OMEGA0*(x[0]*[0,-1]+x[1]*[1,0])
 
78
else:
 
79
   V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.])
 
80
#===================
 
81
fc=TransportPDE(dom,num_equations=1,theta=THETA)
 
82
x=Function(dom).getX()
 
83
fc.setValue(M=Scalar(1.,Function(dom)),C=V)
 
84
#==============
 
85
if TEST_SUPG:
 
86
   supg=LinearSinglePDE(dom)
 
87
   supg.setValue(D=1.)
 
88
   supg.setSolverMethod(supg.LUMPING)
 
89
   dt_supg=inf(dom.getSize()/length(V))
 
90
   u_supg=u0*1.
 
91
 
 
92
c=0
 
93
# saveVTK("u.%s.vtu"%c,u=u0)
 
94
fc.setInitialSolution(u0)
 
95
t=T0
 
96
print "QUALITY FCT: time = %s pi"%(t/pi),inf(u0),sup(u0),integrate(u0)
 
97
while t<T_END:
 
98
    print "time step t=",t+dt   
 
99
    u=fc.solve(dt, verbose=True)
 
100
    print "QUALITY FCT: time = %s pi"%(t+dt/pi),inf(u),sup(u),integrate(u)
 
101
    if TEST_SUPG:
 
102
        #========== supg tests ================
 
103
        nn=max(ceil(dt/dt_supg),1.)
 
104
        dt2=dt/nn
 
105
        nnn=0
 
106
        while nnn<nn :
 
107
            supg.setValue(Y=u_supg+dt2/2*inner(V,grad(u_supg)))
 
108
            u2=supg.getSolution()
 
109
            supg.setValue(Y=u_supg+dt2*inner(V,grad(u2)))
 
110
            u_supg=supg.getSolution()
 
111
            nnn+=1
 
112
    c+=1
 
113
    t+=dt
 
114
    if TEST_SUPG: 
 
115
       print "QUALITY SUPG: time = %s pi"%(t/pi),inf(u_supg),sup(u_supg),integrate(u_supg)
 
116
       # saveVTK("u2.%s.vtu"%c,u=u,u_supg=u_supg)
 
117
    else:
 
118
       # saveVTK("u.%s.vtu"%c,u=u)
 
119
       pass
 
120
    # if c == 20: 1/0