Revision 299 - (show annotations)
Fri Dec 2 05:45:38 2005 UTC (16 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 2841 byte(s)
```runs again
```
 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

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision