/[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 299 - (show annotations)
Fri Dec 2 05:45:38 2005 UTC (13 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 2841 byte(s)
runs again
1 #test_advect.py Gauss
2
3
4 from esys.escript import *
5 import esys.finley
6 from esys.escript.linearPDEs import LinearPDE
7 from esys.escript.linearPDEs import AdvectivePDE
8 import os
9
10 def classical_lhs():
11 ####################### FINLEY UPWINDING HARDCODED #################
12 global surfacePDE
13 surfacePDE = LinearPDE(mesh)
14
15 A = Tensor(0.0, Function(mesh))
16 for j in range(2):
17 for l in range(2):
18 A[j,l] = velocity[j] * velocity[l] * cte * dt
19
20 B = velocity * cte
21 C = dt * velocity
22 D=1.0
23
24 surfacePDE.setValue(A=A, B=B, C=C, D=D)
25
26 def classical_rhs():
27 ####################### FINLEY UPWINDING HARDCODED #################
28 X = velocity * cte * gauss_tronc_old
29 Y = gauss_tronc_old
30 surfacePDE.setValue(X=X, Y=Y)
31
32 def finley_upwd_lhs():
33 ####################### FINLEY UPWINDING #################
34 global surfacePDE
35 surfacePDE = AdvectivePDE(mesh)
36
37 C = dt * velocity
38 D=1.0
39 surfacePDE.setValue(C=C, D=D)
40
41 def finley_upwd_rhs():
42 ####################### FINLEY UPWINDING #################
43 Y = gauss_tronc_old
44 surfacePDE.setValue(Y=Y)
45
46 def taylor_galerkin_lhs():
47 ####################### TAYLOR - GALERKIN #################
48 global surfacePDE
49 surfacePDE = LinearPDE(mesh)
50
51 surfacePDE.setValue(D=1.0)
52
53 def taylor_galerkin_rhs():
54 ####################### TAYLOR - GALERKIN #################
55 X = Vector(0.0, Function(mesh))
56 for i in range(2):
57 for j in range(2):
58 X[i] -= (dt**2)/2.0*velocity[i]*velocity[j]*grad(gauss_tronc_old)[j]
59 Y = gauss_tronc_old*Scalar(1.0, Function(mesh))
60 for j in range(2):
61 Y -= dt*velocity[j]*grad(gauss_tronc_old)[j]
62
63 surfacePDE.setValue(X=X, Y=Y)
64
65
66 ############ start of main code ########################
67
68
69 dt = 0.5
70 t_step = 1
71
72
73 mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)
74
75
76 xx = mesh.getX()[0]
77 yy = mesh.getX()[1]
78
79 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)
80 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
81 gauss_tronc_new = gauss*mask_tronc
82
83 reference = Lsup(gauss_tronc_new)
84
85 h = Lsup(mesh.getSize())
86
87 coeff = 0.3
88
89 v_adim = coeff*h/dt
90 cte = h/(2.0*v_adim)
91 t_step_end = 2.0/(coeff*h)
92
93 velocity = Vector(0.0, Function(mesh))
94 velocity[0] = v_adim
95
96
97 taylor_galerkin_lhs()
98
99 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
100 surfacePDE.setValue(q=q)
101
102
103 while (t_step <= t_step_end):
104
105 gauss_tronc_old = gauss_tronc_new
106
107 taylor_galerkin_rhs()
108
109 gauss_tronc_new = surfacePDE.getSolution()
110
111 print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
112
113 dist = v_adim*dt*t_step
114
115 error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
116 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