/[escript]/temp_trunk_copy/modellib/py_src/flow.py
ViewVC logotype

Diff of /temp_trunk_copy/modellib/py_src/flow.py

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

trunk/esys2/modellib/py_src/flow.py revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC trunk/modellib/py_src/flow.py revision 918 by gross, Wed Jan 3 06:30:00 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3  from escript.escript import *  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
4  from escript.modelframe import Model,IterationDivergenceError                      http://www.access.edu.au
5  from escript.linearPDEs import LameEquation                  Primary Business: Queensland, Australia"""
6    __license__="""Licensed under the Open Software License version 3.0
7                 http://www.opensource.org/licenses/osl-3.0.php"""
8    
9    from esys.escript import *
10    from esys.escript.modelframe import Model,IterationDivergenceError
11    from esys.escript.linearPDEs import LameEquation
12    
13  class SteadyIncompressibleFlow(Model):  class SteadyIncompressibleFlow(Model):
14         """         """
15    
16                     \f[         M{-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i}
17                     -\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i        
18                     \f]         M{\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}}
19    
20                     \f[         M{D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right)}
21                     \sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}        
22                     \f[         M{v_{k,k} = 0}
                    D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; .  
                    \f]  
   
                    \f[  
                    -v_{k,k} = 0 \; .  
                    \f]  
23    
24         """         """
25    
26         def __init__(self,debug=False):         def __init__(self,**kwargs):
27             """set up model"""             """
28             Model.__init__(self,debug=debug)             set up model
29               """
30               Model.__init__(self,**kwargs)
31             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
32                                   velocity=0., \                                   velocity=0., \
33                                   pressure=0., \                                   pressure=0., \
34                                   viscosity=1., \                                   viscosity=1., \
35                                   internal_force=0., \                                   internal_force=0., \
36                                   location_prescribed_velocity=Data(), \                                   location_prescribed_velocity=None, \
37                                   prescribed_velocity=Data(), \                                   prescribed_velocity=None, \
38                                   rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.001)                                   rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.0001)
39             self.__iter=0             self.__iter=0
40    
41         def doInitialization(self):         def doInitialization(self):
42             """initialize model"""             """
43               initialize model
44               """
45             self.__p_old=None             self.__p_old=None
46             self.__p_very_old=None             self.__p_very_old=None
47             self.__dt_old=None             self.__dt_old=None
48             self.__pde=LameEquation(self.domain)             self.__pde=LameEquation(self.domain)
49               self.__pde.setSolverMethod(LameEquation.DIRECT)
50               if self.location_prescribed_velocity == None: self.location_prescribed_velocit=Data()
51               if self.prescribed_velocity == None: self.prescribed_velocity=Data()
52    
53         def stress(self):         def stress(self):
54             """returns current stress"""             """
55               returns current stress"""
56             return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain)             return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain)
57    
58         def stretching(self):         def stretching(self):
59             """returns stertching tensor"""             """
60               returns stertching tensor
61               """
62             g=grad(self.velocity)             g=grad(self.velocity)
63             return (g+transpose(g))/2             return (g+transpose(g))/2
64    
65         def doStepPreprocessing(self,dt):         def doStepPreprocessing(self,dt):
66              """              """
67                 step up pressure iteration              step up pressure iteration
68    
69                 if run within a time dependend problem extrapolation of pressure from previous time steps is used to              if run within a time dependend problem extrapolation of pressure from previous time steps is used to
70                 get an initial guess (that needs some work!!!!!!!)              get an initial guess (that needs some work!!!!!!!)
71              """              """
72              self.__iter=0              self.__iter=0
73              self.__diff=1.e40              self.__diff=1.e40
# Line 70  class SteadyIncompressibleFlow(Model): Line 80  class SteadyIncompressibleFlow(Model):
80         def doStep(self,dt):         def doStep(self,dt):
81            """            """
82    
83               performs an iteration step of the penalty method.            performs an iteration step of the penalty method.
84               IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.            IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
85    
86            """            """
87            penalty=self.viscosity/self.relaxation            penalty=self.viscosity/self.relaxation
# Line 84  class SteadyIncompressibleFlow(Model): Line 94  class SteadyIncompressibleFlow(Model):
94            self.__pde.setTolerance(self.rel_tol/10.)            self.__pde.setTolerance(self.rel_tol/10.)
95            self.velocity=self.__pde.getSolution()            self.velocity=self.__pde.getSolution()
96            update=penalty*div(self.velocity)            update=penalty*div(self.velocity)
           print "f=",inf(update),sup(update)  
97            self.pressure=self.pressure-update            self.pressure=self.pressure-update
98            self.__diff,diff_old=Lsup(update),self.__diff            self.__diff,diff_old=Lsup(update),self.__diff
99            self.__iter+=1            self.__iter+=1

Legend:
Removed from v.147  
changed lines
  Added in v.918

  ViewVC Help
Powered by ViewVC 1.1.26