/[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 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 3241 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 #test_advect.py Gauss
2
3
4 from esys.escript import *
5 import esys.finley
6 from esys.linearPDEs import LinearPDE
7 from esys.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 if os.path.exists("results/sylvain/dist_error.dat"):
76 os.unlink("results/sylvain/dist_error.dat")
77
78
79 xx = mesh.getX()[0]
80 yy = mesh.getX()[1]
81
82 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)
83 mask_tronc = (xx-0.25).wherePositive()*(1.25-xx).wherePositive()*(yy-0.5).wherePositive()*(1.5-yy).wherePositive()
84 gauss_tronc_new = gauss*mask_tronc
85
86 reference = Lsup(gauss_tronc_new)
87 gauss_tronc_new.saveDX("results/sylvain/gauss.%2.2i.dx" % 0)
88
89 h = Lsup(mesh.getSize())
90
91 coeff = 0.3
92
93 v_adim = coeff*h/dt
94 cte = h/(2.0*v_adim)
95 t_step_end = 2.0/(coeff*h)
96
97 velocity = Vector(0.0, Function(mesh))
98 velocity[0] = v_adim
99 velocity.saveDX("results/sylvain/velocity_field.dx")
100
101
102 taylor_galerkin_lhs()
103
104 q = mesh.getX()[0].whereZero() + (4.0-mesh.getX()[0]).whereZero() + mesh.getX()[1].whereZero() + (2.0-mesh.getX()[1]).whereZero()
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 gauss_tronc_old.saveDX("results/sylvain/gauss.%2.2i.dx" % t_step)
117 print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
118
119 dist = v_adim*dt*t_step
120
121 error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
122 File1 = file("results/sylvain/dist_error.dat", "a", 1)
123 File1.write("%e %e\n" % (dist, error))
124 File1.close()
125 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