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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 299 - (hide annotations)
Fri Dec 2 05:45:38 2005 UTC (13 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 2841 byte(s)
runs again
1 jgs 123 #test_advect.py Gauss
2    
3    
4     from esys.escript import *
5     import esys.finley
6 gross 299 from esys.escript.linearPDEs import LinearPDE
7     from esys.escript.linearPDEs import AdvectivePDE
8 jgs 123 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 gross 299 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
81 jgs 123 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 gross 299 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
100 jgs 123 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