| 1 |
#test_advect.py Gauss |
| 2 |
|
| 3 |
|
| 4 |
from esys.escript import * |
| 5 |
import esys.finley |
| 6 |
from esys.escript.linearPDEs import LinearPDE |
| 7 |
from esys.escript.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.0*velocity[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 -= dt*velocity[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 |
|
| 76 |
xx = mesh.getX()[0] |
| 77 |
yy = mesh.getX()[1] |
| 78 |
|
| 79 |
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) |
| 80 |
mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy) |
| 81 |
gauss_tronc_new = gauss*mask_tronc |
| 82 |
|
| 83 |
reference = Lsup(gauss_tronc_new) |
| 84 |
|
| 85 |
h = Lsup(mesh.getSize()) |
| 86 |
|
| 87 |
coeff = 0.3 |
| 88 |
|
| 89 |
v_adim = coeff*h/dt |
| 90 |
cte = h/(2.0*v_adim) |
| 91 |
t_step_end = 2.0/(coeff*h) |
| 92 |
|
| 93 |
velocity = Vector(0.0, Function(mesh)) |
| 94 |
velocity[0] = v_adim |
| 95 |
|
| 96 |
|
| 97 |
taylor_galerkin_lhs() |
| 98 |
|
| 99 |
q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1]) |
| 100 |
surfacePDE.setValue(q=q) |
| 101 |
|
| 102 |
|
| 103 |
while (t_step <= t_step_end): |
| 104 |
|
| 105 |
gauss_tronc_old = gauss_tronc_new |
| 106 |
|
| 107 |
taylor_galerkin_rhs() |
| 108 |
|
| 109 |
gauss_tronc_new = surfacePDE.getSolution() |
| 110 |
|
| 111 |
print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh))) |
| 112 |
|
| 113 |
dist = v_adim*dt*t_step |
| 114 |
|
| 115 |
error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference |
| 116 |
t_step = t_step + 1 |
Parent Directory
| 
Revision Log