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