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

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

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

trunk/esys2/finley/test/python/test_advect.py revision 123 by jgs, Fri Jul 8 04:08:13 2005 UTC trunk/finley/test/python/test_advect.py revision 617 by elspeth, Wed Mar 22 02:58:17 2006 UTC
# Line 1  Line 1 
1  #test_advect.py Gauss  #test_advect.py Gauss
2    
3    
4    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
5                        http://www.access.edu.au
6                    Primary Business: Queensland, Australia"""
7    __license__="""Licensed under the Open Software License version 3.0
8                 http://www.opensource.org/licenses/osl-3.0.php"""
9  from esys.escript import *  from esys.escript import *
10  import esys.finley  import esys.finley
11  from esys.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
12  from esys.linearPDEs import AdvectivePDE  from esys.escript.linearPDEs import AdvectivePDE
13  import os  import os
14    
15  def classical_lhs():  def classical_lhs():
# Line 72  t_step = 1 Line 77  t_step = 1
77    
78  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)
79    
 if os.path.exists("results/sylvain/dist_error.dat"):  
   os.unlink("results/sylvain/dist_error.dat")  
   
80    
81  xx = mesh.getX()[0]  xx = mesh.getX()[0]
82  yy = mesh.getX()[1]  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)  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 = (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)
86  gauss_tronc_new = gauss*mask_tronc  gauss_tronc_new = gauss*mask_tronc
87    
88  reference = Lsup(gauss_tronc_new)  reference = Lsup(gauss_tronc_new)
 gauss_tronc_new.saveDX("results/sylvain/gauss.%2.2i.dx" % 0)  
89    
90  h = Lsup(mesh.getSize())  h = Lsup(mesh.getSize())
91    
# Line 96  t_step_end = 2.0/(coeff*h) Line 97  t_step_end = 2.0/(coeff*h)
97    
98  velocity = Vector(0.0, Function(mesh))  velocity = Vector(0.0, Function(mesh))
99  velocity[0] = v_adim  velocity[0] = v_adim
 velocity.saveDX("results/sylvain/velocity_field.dx")  
100    
101    
102  taylor_galerkin_lhs()  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()  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)  surfacePDE.setValue(q=q)
106    
107    
# Line 113  while (t_step <= t_step_end): Line 113  while (t_step <= t_step_end):
113        
114    gauss_tronc_new = surfacePDE.getSolution()    gauss_tronc_new = surfacePDE.getSolution()
115    
   gauss_tronc_old.saveDX("results/sylvain/gauss.%2.2i.dx" % t_step)  
116    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)))
117        
118    dist = v_adim*dt*t_step    dist = v_adim*dt*t_step
119        
120    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()  
121    t_step = t_step + 1    t_step = t_step + 1

Legend:
Removed from v.123  
changed lines
  Added in v.617

  ViewVC Help
Powered by ViewVC 1.1.26