/[escript]/trunk/finley/test/python/test_advect.py
ViewVC logotype

Diff of /trunk/finley/test/python/test_advect.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 298 by jgs, Wed Nov 9 02:02:19 2005 UTC revision 299 by gross, Fri Dec 2 05:45:38 2005 UTC
# Line 3  Line 3 
3    
4  from esys.escript import *  from esys.escript import *
5  import esys.finley  import esys.finley
6  from esys.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
7  from esys.linearPDEs import AdvectivePDE  from esys.escript.linearPDEs import AdvectivePDE
8  import os  import os
9    
10  def classical_lhs():  def classical_lhs():
# Line 72  t_step = 1 Line 72  t_step = 1
72    
73  mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)  mesh = esys.finley.Rectangle(l0=4.0, l1=2.0, order=1, n0=60, n1=30)
74    
 if os.path.exists("results/sylvain/dist_error.dat"):  
   os.unlink("results/sylvain/dist_error.dat")  
   
75    
76  xx = mesh.getX()[0]  xx = mesh.getX()[0]
77  yy = mesh.getX()[1]  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)  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 = (xx-0.25).wherePositive()*(1.25-xx).wherePositive()*(yy-0.5).wherePositive()*(1.5-yy).wherePositive()  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  gauss_tronc_new = gauss*mask_tronc
82    
83  reference = Lsup(gauss_tronc_new)  reference = Lsup(gauss_tronc_new)
 gauss_tronc_new.saveDX("results/sylvain/gauss.%2.2i.dx" % 0)  
84    
85  h = Lsup(mesh.getSize())  h = Lsup(mesh.getSize())
86    
# Line 96  t_step_end = 2.0/(coeff*h) Line 92  t_step_end = 2.0/(coeff*h)
92    
93  velocity = Vector(0.0, Function(mesh))  velocity = Vector(0.0, Function(mesh))
94  velocity[0] = v_adim  velocity[0] = v_adim
 velocity.saveDX("results/sylvain/velocity_field.dx")  
95    
96    
97  taylor_galerkin_lhs()  taylor_galerkin_lhs()
98    
99  q = mesh.getX()[0].whereZero() + (4.0-mesh.getX()[0]).whereZero() + mesh.getX()[1].whereZero() + (2.0-mesh.getX()[1]).whereZero()  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)  surfacePDE.setValue(q=q)
101    
102    
# Line 113  while (t_step <= t_step_end): Line 108  while (t_step <= t_step_end):
108        
109    gauss_tronc_new = surfacePDE.getSolution()    gauss_tronc_new = surfacePDE.getSolution()
110    
   gauss_tronc_old.saveDX("results/sylvain/gauss.%2.2i.dx" % t_step)  
111    print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))    print "integral of f", integrate(gauss_tronc_new*Scalar(1.0, Function(mesh)))
112        
113    dist = v_adim*dt*t_step    dist = v_adim*dt*t_step
114        
115    error = 100*abs(Lsup(gauss_tronc_new)-reference)/reference    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()  
116    t_step = t_step + 1    t_step = t_step + 1
Legend:
Removed from v.298  
changed lines
  Added in v.299
  ViewVC Help
Powered by ViewVC 1.1.26