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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/finley/test/python/test_advect.py
File MIME type: text/x-python
File size: 3241 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 jgs 123 #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