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

Contents of /trunk/finley/test/python/test_advect.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/x-python
File size: 3129 byte(s)
More copyright.

1 #test_advect.py Gauss
2
3
4 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
5 http://www.access.edu.au
6 Primary Business: Queensland, Australia"""
7 __license__="""Licensed under the Open Software License version 3.0
8 http://www.opensource.org/licenses/osl-3.0.php"""
9 from esys.escript import *
10 import esys.finley
11 from esys.escript.linearPDEs import LinearPDE
12 from esys.escript.linearPDEs import AdvectivePDE
13 import os
14
15 def classical_lhs():
16 ####################### FINLEY UPWINDING HARDCODED #################
17 global surfacePDE
18 surfacePDE = LinearPDE(mesh)
19
20 A = Tensor(0.0, Function(mesh))
21 for j in range(2):
22 for l in range(2):
23 A[j,l] = velocity[j] * velocity[l] * cte * dt
24
25 B = velocity * cte
26 C = dt * velocity
27 D=1.0
28
29 surfacePDE.setValue(A=A, B=B, C=C, D=D)
30
31 def classical_rhs():
32 ####################### FINLEY UPWINDING HARDCODED #################
33 X = velocity * cte * gauss_tronc_old
34 Y = gauss_tronc_old
35 surfacePDE.setValue(X=X, Y=Y)
36
37 def finley_upwd_lhs():
38 ####################### FINLEY UPWINDING #################
39 global surfacePDE
40 surfacePDE = AdvectivePDE(mesh)
41
42 C = dt * velocity
43 D=1.0
44 surfacePDE.setValue(C=C, D=D)
45
46 def finley_upwd_rhs():
47 ####################### FINLEY UPWINDING #################
48 Y = gauss_tronc_old
49 surfacePDE.setValue(Y=Y)
50
51 def taylor_galerkin_lhs():
52 ####################### TAYLOR - GALERKIN #################
53 global surfacePDE
54 surfacePDE = LinearPDE(mesh)
55
56 surfacePDE.setValue(D=1.0)
57
58 def taylor_galerkin_rhs():
59 ####################### TAYLOR - GALERKIN #################
60 X = Vector(0.0, Function(mesh))
61 for i in range(2):
62 for j in range(2):
63 X[i] -= (dt**2)/2.0*velocity[i]*velocity[j]*grad(gauss_tronc_old)[j]
64 Y = gauss_tronc_old*Scalar(1.0, Function(mesh))
65 for j in range(2):
66 Y -= dt*velocity[j]*grad(gauss_tronc_old)[j]
67
68 surfacePDE.setValue(X=X, Y=Y)
69
70
71 ############ start of main code ########################
72
73
74 dt = 0.5
75 t_step = 1
76
77
78 mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)
79
80
81 xx = mesh.getX()[0]
82 yy = mesh.getX()[1]
83
84 gauss = (160*(xx-0.25)**4 - 320*(xx-0.25)**3 + 160*(xx-0.25)**2 )*( 160*(yy-0.5)**4 - 320*(yy-0.5)**3 + 160*(yy-0.5)**2)
85 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
86 gauss_tronc_new = gauss*mask_tronc
87
88 reference = Lsup(gauss_tronc_new)
89
90 h = Lsup(mesh.getSize())
91
92 coeff = 0.3
93
94 v_adim = coeff*h/dt
95 cte = h/(2.0*v_adim)
96 t_step_end = 2.0/(coeff*h)
97
98 velocity = Vector(0.0, Function(mesh))
99 velocity[0] = v_adim
100
101
102 taylor_galerkin_lhs()
103
104 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
105 surfacePDE.setValue(q=q)
106
107
108 while (t_step <= t_step_end):
109
110 gauss_tronc_old = gauss_tronc_new
111
112 taylor_galerkin_rhs()
113
114 gauss_tronc_new = surfacePDE.getSolution()
115
116 print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
117
118 dist = v_adim*dt*t_step
119
120 error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
121 t_step = t_step + 1

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26