test/python/advective.py

#
# $Id$
#
#######################################################
#
#           Copyright 2003-2007 by ACceSS MNRF
#       Copyright 2007 by University of Queensland
#
#                http://esscc.uq.edu.au
#        Primary Business: Queensland, Australia
#  Licensed under the Open Software License version 3.0
#     http://www.opensource.org/licenses/osl-3.0.php
#
#######################################################
#

#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	#
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.0velocity[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 -= dtvelocity[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_adimdtt_step
135
136	error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
137	t_step = t_step + 1
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision