/[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 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 2 months ago) by elspeth
Original Path: trunk/finley/test/python/test_advect.py
File MIME type: text/x-python
File size: 3129 byte(s)
More copyright.

1 jgs 123 #test_advect.py Gauss
2    
3    
4 elspeth 617 __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 jgs 123 from esys.escript import *
10     import esys.finley
11 gross 299 from esys.escript.linearPDEs import LinearPDE
12     from esys.escript.linearPDEs import AdvectivePDE
13 jgs 123 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 gross 299 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
86 jgs 123 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 gross 299 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
105 jgs 123 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