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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2549 - (show annotations)
Mon Jul 20 06:43:47 2009 UTC (10 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 3600 byte(s)
Remainder of copyright date fixes
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 #test_advect.py Gauss
23
24 from esys.escript import *
25 import esys.finley
26 from esys.escript.linearPDEs import LinearPDE
27 from esys.escript.linearPDEs import AdvectivePDE
28 import os
29
30 def classical_lhs():
31 ####################### FINLEY UPWINDING HARDCODED #################
32 global surfacePDE
33 surfacePDE = LinearPDE(mesh)
34
35 A = Tensor(0.0, Function(mesh))
36 for j in range(2):
37 for l in range(2):
38 A[j,l] = velocity[j] * velocity[l] * cte * dt
39
40 B = velocity * cte
41 C = dt * velocity
42 D=1.0
43
44 surfacePDE.setValue(A=A, B=B, C=C, D=D)
45
46 def classical_rhs():
47 ####################### FINLEY UPWINDING HARDCODED #################
48 X = velocity * cte * gauss_tronc_old
49 Y = gauss_tronc_old
50 surfacePDE.setValue(X=X, Y=Y)
51
52 def finley_upwd_lhs():
53 ####################### FINLEY UPWINDING #################
54 global surfacePDE
55 surfacePDE = AdvectivePDE(mesh)
56
57 C = dt * velocity
58 D=1.0
59 surfacePDE.setValue(C=C, D=D)
60
61 def finley_upwd_rhs():
62 ####################### FINLEY UPWINDING #################
63 Y = gauss_tronc_old
64 surfacePDE.setValue(Y=Y)
65
66 def taylor_galerkin_lhs():
67 ####################### TAYLOR - GALERKIN #################
68 global surfacePDE
69 surfacePDE = LinearPDE(mesh)
70
71 surfacePDE.setValue(D=1.0)
72
73 def taylor_galerkin_rhs():
74 ####################### TAYLOR - GALERKIN #################
75 X = Vector(0.0, Function(mesh))
76 for i in range(2):
77 for j in range(2):
78 X[i] -= (dt**2)/2.0*velocity[i]*velocity[j]*grad(gauss_tronc_old)[j]
79 Y = gauss_tronc_old*Scalar(1.0, Function(mesh))
80 for j in range(2):
81 Y -= dt*velocity[j]*grad(gauss_tronc_old)[j]
82
83 surfacePDE.setValue(X=X, Y=Y)
84
85
86 ############ start of main code ########################
87
88
89 dt = 0.5
90 t_step = 1
91
92
93 mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)
94
95
96 xx = mesh.getX()[0]
97 yy = mesh.getX()[1]
98
99 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)
100 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
101 gauss_tronc_new = gauss*mask_tronc
102
103 reference = Lsup(gauss_tronc_new)
104
105 h = Lsup(mesh.getSize())
106
107 coeff = 0.3
108
109 v_adim = coeff*h/dt
110 cte = h/(2.0*v_adim)
111 t_step_end = 2.0/(coeff*h)
112
113 velocity = Vector(0.0, Function(mesh))
114 velocity[0] = v_adim
115
116
117 taylor_galerkin_lhs()
118
119 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
120 surfacePDE.setValue(q=q)
121
122
123 while (t_step <= t_step_end):
124
125 gauss_tronc_old = gauss_tronc_new
126
127 taylor_galerkin_rhs()
128
129 gauss_tronc_new = surfacePDE.getSolution()
130
131 print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
132
133 dist = v_adim*dt*t_step
134
135 error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
136 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