Revision 617 - (show annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years, 2 months ago) by elspeth
File MIME type: text/x-python
File size: 3129 byte(s)
```More copyright.

```
 1 #test_advect.py Gauss 2 3 4 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF 5 6 Primary Business: Queensland, Australia""" 7 __license__="""Licensed under the Open Software License version 3.0 8 9 from esys.escript import * 10 import esys.finley 11 from esys.escript.linearPDEs import LinearPDE 12 from esys.escript.linearPDEs import AdvectivePDE 13 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.0*velocity[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 -= dt*velocity[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 mask_tronc = wherePositive(xx-0.25)*wherePositive(1.25-xx)*wherePositive(yy-0.5)*wherePositive(1.5-yy) 86 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 q = whereZero(mesh.getX()[0]) + whereZero(4.0-mesh.getX()[0]) + whereZero(mesh.getX()[1]) + whereZero(2.0-mesh.getX()[1]) 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 print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh))) 117 118 dist = v_adim*dt*t_step 119 120 error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference 121 t_step = t_step + 1

## Properties

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