test/python/advective.py

#test_advect.py Gauss


__copyright__="""  Copyright (c) 2006 by ACcESS MNRF
                    http://www.access.edu.au
                Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
             http://www.opensource.org/licenses/osl-3.0.php"""
from esys.escript import *
import esys.finley
from esys.escript.linearPDEs import LinearPDE
from esys.escript.linearPDEs import AdvectivePDE
import os

def classical_lhs():
#######################   FINLEY UPWINDING HARDCODED    #################
  global surfacePDE
  surfacePDE = LinearPDE(mesh)
  
  A = Tensor(0.0, Function(mesh))
  for j in range(2):
    for l in range(2):
      A[j,l] = velocity[j] * velocity[l] * cte * dt
     
  B = velocity * cte
  C = dt * velocity
  D=1.0
  
  surfacePDE.setValue(A=A, B=B, C=C,  D=D)

def classical_rhs():
#######################   FINLEY UPWINDING HARDCODED    #################
  X = velocity * cte * gauss_tronc_old
  Y = gauss_tronc_old
  surfacePDE.setValue(X=X, Y=Y)
  
def finley_upwd_lhs():
#######################   FINLEY UPWINDING   #################
  global surfacePDE
  surfacePDE = AdvectivePDE(mesh)
  
  C = dt * velocity
  D=1.0
  surfacePDE.setValue(C=C,  D=D)
  
def finley_upwd_rhs():
#######################   FINLEY UPWINDING   #################
  Y = gauss_tronc_old
  surfacePDE.setValue(Y=Y)

def taylor_galerkin_lhs():
#######################   TAYLOR - GALERKIN   #################
  global surfacePDE
  surfacePDE = LinearPDE(mesh)
  
  surfacePDE.setValue(D=1.0)
  
def taylor_galerkin_rhs():
#######################   TAYLOR - GALERKIN   #################
  X = Vector(0.0, Function(mesh))
  for i in range(2):
    for j in range(2):
      X[i] -= (dt**2)/2.0*velocity[i]*velocity[j]*grad(gauss_tronc_old)[j]
  Y = gauss_tronc_old*Scalar(1.0, Function(mesh))
  for j in range(2):
    Y -= dt*velocity[j]*grad(gauss_tronc_old)[j]
    
  surfacePDE.setValue(X=X, Y=Y)


############ start of main code ########################


dt = 0.5
t_step = 1


mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)


xx = mesh.getX()[0]
yy = mesh.getX()[1]

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)
mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy)
gauss_tronc_new = gauss*mask_tronc

reference = Lsup(gauss_tronc_new)

h = Lsup(mesh.getSize())

coeff = 0.3

v_adim = coeff*h/dt
cte = h/(2.0*v_adim)
t_step_end = 2.0/(coeff*h)

velocity = Vector(0.0, Function(mesh))
velocity[0] = v_adim


taylor_galerkin_lhs()

q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1])
surfacePDE.setValue(q=q)


while (t_step <= t_step_end):
   
  gauss_tronc_old = gauss_tronc_new
  
  taylor_galerkin_rhs()
  
  gauss_tronc_new = surfacePDE.getSolution()

  print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
  
  dist = v_adim*dt*t_step
  
  error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
  t_step = t_step + 1
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.0velocity[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 -= dtvelocity[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_adimdtt_step
119
120			error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
121			t_step = t_step + 1
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision