test/python/test_advect.py

#test_advect.py Gauss


from esys.escript import *
import esys.finley
from esys.linearPDEs import LinearPDE
from esys.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)

if os.path.exists("results/sylvain/dist_error.dat"):
  os.unlink("results/sylvain/dist_error.dat")


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

reference = Lsup(gauss_tronc_new)
gauss_tronc_new.saveDX("results/sylvain/gauss.%2.2i.dx" % 0)

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
velocity.saveDX("results/sylvain/velocity_field.dx")


taylor_galerkin_lhs()

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


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

  gauss_tronc_old.saveDX("results/sylvain/gauss.%2.2i.dx" % t_step)
  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
  File1 = file("results/sylvain/dist_error.dat", "a", 1)
  File1.write("%e %e\n" % (dist, error))
  File1.close()
  t_step = t_step + 1
1	jgs	123	#test_advect.py Gauss
2
3
4			from esys.escript import *
5			import esys.finley
6			from esys.linearPDEs import LinearPDE
7			from esys.linearPDEs import AdvectivePDE
8			import os
9
10			def classical_lhs():
11			####################### FINLEY UPWINDING HARDCODED #################
12			global surfacePDE
13			surfacePDE = LinearPDE(mesh)
14
15			A = Tensor(0.0, Function(mesh))
16			for j in range(2):
17			for l in range(2):
18			A[j,l] = velocity[j] * velocity[l] * cte * dt
19
20			B = velocity * cte
21			C = dt * velocity
22			D=1.0
23
24			surfacePDE.setValue(A=A, B=B, C=C, D=D)
25
26			def classical_rhs():
27			####################### FINLEY UPWINDING HARDCODED #################
28			X = velocity * cte * gauss_tronc_old
29			Y = gauss_tronc_old
30			surfacePDE.setValue(X=X, Y=Y)
31
32			def finley_upwd_lhs():
33			####################### FINLEY UPWINDING #################
34			global surfacePDE
35			surfacePDE = AdvectivePDE(mesh)
36
37			C = dt * velocity
38			D=1.0
39			surfacePDE.setValue(C=C, D=D)
40
41			def finley_upwd_rhs():
42			####################### FINLEY UPWINDING #################
43			Y = gauss_tronc_old
44			surfacePDE.setValue(Y=Y)
45
46			def taylor_galerkin_lhs():
47			####################### TAYLOR - GALERKIN #################
48			global surfacePDE
49			surfacePDE = LinearPDE(mesh)
50
51			surfacePDE.setValue(D=1.0)
52
53			def taylor_galerkin_rhs():
54			####################### TAYLOR - GALERKIN #################
55			X = Vector(0.0, Function(mesh))
56			for i in range(2):
57			for j in range(2):
58			X[i] -= (dt*2)/2.0velocity[i]velocity[j]grad(gauss_tronc_old)[j]
59			Y = gauss_tronc_old*Scalar(1.0, Function(mesh))
60			for j in range(2):
61			Y -= dtvelocity[j]grad(gauss_tronc_old)[j]
62
63			surfacePDE.setValue(X=X, Y=Y)
64
65
66			############ start of main code ########################
67
68
69			dt = 0.5
70			t_step = 1
71
72
73			mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)
74
75			if os.path.exists("results/sylvain/dist_error.dat"):
76			os.unlink("results/sylvain/dist_error.dat")
77
78
79			xx = mesh.getX()[0]
80			yy = mesh.getX()[1]
81
82			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)
83			mask_tronc = (xx-0.25).wherePositive()(1.25-xx).wherePositive()(yy-0.5).wherePositive()*(1.5-yy).wherePositive()
84			gauss_tronc_new = gauss*mask_tronc
85
86			reference = Lsup(gauss_tronc_new)
87			gauss_tronc_new.saveDX("results/sylvain/gauss.%2.2i.dx" % 0)
88
89			h = Lsup(mesh.getSize())
90
91			coeff = 0.3
92
93			v_adim = coeff*h/dt
94			cte = h/(2.0*v_adim)
95			t_step_end = 2.0/(coeff*h)
96
97			velocity = Vector(0.0, Function(mesh))
98			velocity[0] = v_adim
99			velocity.saveDX("results/sylvain/velocity_field.dx")
100
101
102			taylor_galerkin_lhs()
103
104			q = mesh.getX()[0].whereZero() + (4.0-mesh.getX()[0]).whereZero() + mesh.getX()[1].whereZero() + (2.0-mesh.getX()[1]).whereZero()
105			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			gauss_tronc_old.saveDX("results/sylvain/gauss.%2.2i.dx" % t_step)
117			print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
118
119			dist = v_adimdtt_step
120
121			error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference
122			File1 = file("results/sylvain/dist_error.dat", "a", 1)
123			File1.write("%e %e\n" % (dist, error))
124			File1.close()
125			t_step = t_step + 1
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision