/[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 144 by jgs, Mon Jul 25 05:42:21 2005 UTC trunk/modellib/py_src/flow.py revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
11    #  Licensed under the Open Software License version 3.0
12    #     http://www.opensource.org/licenses/osl-3.0.php
13    #
14    #######################################################
15    #
16    
17    
18    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
19                        http://www.access.edu.au
20                    Primary Business: Queensland, Australia"""
21    __license__="""Licensed under the Open Software License version 3.0
22                 http://www.opensource.org/licenses/osl-3.0.php"""
23    
24  from esys.escript import *  from esys.escript import *
25  from esys.modelframe import Model,IterationDivergenceError  from esys.escript.modelframe import Model,IterationDivergenceError
26  from esys.linearPDEs import LameEquation  from esys.escript.linearPDEs import LameEquation
27    
28  class SteadyIncompressibleFlow(Model):  class SteadyIncompressibleFlow(Model):
29         """         """
30    
31           M{-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i}
32          
33           M{\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}}
34    
35           M{D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right)}
36          
37           M{v_{k,k} = 0}
38    
                    \f[  
                    -\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i  
                    \f]  
   
                    \f[  
                    \sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}  
                    \f[  
                    D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; .  
                    \f]  
   
                    \f[  
                    -v_{k,k} = 0 \; .  
                    \f]  
                     
39         """         """
40    
41         def __init__(self,debug=False):         def __init__(self,**kwargs):
42             """set up model"""             """
43             Model.__init__(self,debug=debug)             set up model
44               """
45               Model.__init__(self,**kwargs)
46             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
47                                   velocity=0., \                                   velocity=0., \
48                                   pressure=0., \                                   pressure=0., \
49                                   viscocity=1., \                                   viscosity=1., \
50                                   internal_force=0., \                                   internal_force=0., \
51                                   location_prescribed_velocities=Data(), \                                   location_prescribed_velocity=None, \
52                                   prescribed_velocities=Data(), \                                   prescribed_velocity=None, \
53                                   rel_tol=1.e-8,abs_tol=1.e-15,max_iter=10,relaxation=0.01)                                   rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.0001)
54             self.__iter=0             self.__iter=0
55    
56         def doInitialization(self,t):         def doInitialization(self):
57             """initialize model"""             """
58             self.__pressure_history=None             initialize model
59             self.__dt_history=None             """
60               self.__p_old=None
61               self.__p_very_old=None
62               self.__dt_old=None
63             self.__pde=LameEquation(self.domain)             self.__pde=LameEquation(self.domain)
64               self.__pde.setSolverMethod(LameEquation.DIRECT)
65               if self.location_prescribed_velocity == None: self.location_prescribed_velocit=Data()
66               if self.prescribed_velocity == None: self.prescribed_velocity=Data()
67    
68         def stress(self):         def stress(self):
69             """returns current stress"""             """
70             return 2*self.viscocity*self.stretching-self.pressure             returns current stress"""
71               return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain)
72    
73         def stretching(self):         def stretching(self):
74             """returns stertching tensor"""             """
75               returns stertching tensor
76               """
77             g=grad(self.velocity)             g=grad(self.velocity)
78             return (g+transpose(g))/2             return (g+transpose(g))/2
79              
80         def doIterationInitialization(self,dt):         def doStepPreprocessing(self,dt):
81              """              """
82                 step up pressure iteration              step up pressure iteration
83    
84                 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
85                 get an initial guess              get an initial guess (that needs some work!!!!!!!)
86              """              """
87              if self.__dt_history==None:              self.__iter=0
88                 self.__pressure_history=self.pressure              self.__diff=1.e40
89                 self.__iter=0              if not self.__p_old==None:
90              else:                 if self.__p_very_old==None:
91                 self.__pressure_history,self.pressure=self.pressure,(1.+dt/self.__dt_history)*self.pressure+dt/self.__dt_history*self.__pressure_history                    self.pressure=self.__p_old
92                 self.__iter=1                 else:
93              self.diff=1.e400                    self.pressure=(1.+dt/self.__dt_old)*self.__p_old-dt/self.__dt_old*self.__p_very_old
94    
95         def doIterationStep(self,dt):         def doStep(self,dt):
96            """            """
97    
98               performs an iteration step of the penalty method.            performs an iteration step of the penalty method.
99               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.
100    
101            """            """
102            penalty=self.viscosity/self.relaxation            penalty=self.viscosity/self.relaxation
103            self.__pde.setValue(lame_lambda=self.viscosity, \            self.__pde.setValue(lame_mu=self.viscosity, \
104                                lame_my=penalty, \                                lame_lambda=penalty, \
105                                F=F, \                                F=self.internal_force, \
106                                sigma0=self.pressure*kronecker(self.__pde.getDomain()), \                                sigma=self.pressure*kronecker(self.__pde.getDomain()), \
107                                r=self.location_prescribed_velocities, \                                r=self.prescribed_velocity, \
108                                q=self.location_prescribed_velocities)                                q=self.location_prescribed_velocity)
109            self.__pde.setTolerance(self.tol*1.e-2)            self.__pde.setTolerance(self.rel_tol/10.)
110            self.velocity=self.pde.getSolution()            self.velocity=self.__pde.getSolution()
111            update=penalty*div(self.velocity)            update=penalty*div(self.velocity)
112            self.pressure=self.pressure-update            self.pressure=self.pressure-update
113            self.diff,diff=Lsup(update),self.diff            self.__diff,diff_old=Lsup(update),self.__diff
114            print "Pressure iteration: step %d: correction %e"%(self.__iter,self.diff/Lsup(self.pressure))            self.__iter+=1
115            if diff<=self.diff:            self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
116              self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
117              self.trace("pressure correction: %e"%self.__diff)
118              if self.__iter>2 and diff_old<self.__diff:
119                  self.trace("Pressure iteration failed!")
120                raise IterationDivergenceError,"no improvement in pressure iteration"                raise IterationDivergenceError,"no improvement in pressure iteration"
121              if self.__iter>self.max_iter:
122                  raise IterationDivergenceError,"Maximum number of iterations steps reached"
123    
124         def terminate(self):         def terminateIteration(self):
125            """iteration is terminated if relative pressure change is less then rel_tol"""            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
126            if self.iter<2:            return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol
127                return False  
128            else:         def doStepPostprocessing(self,dt):
129               return self.diff<self.rel_tol*Lsup(self.pressure)+self.abs_tol            self.__dt_old=dt
130              self.__p_old,self.__p_very_old=self.pressure,self.__p_old
        def doIterationFinalization(self,dt):  
           self.__dt_history=dt  
 if __name__=="__main__":  
   
   
   
    from esys.modelframe import Link,Simulation,ExplicitSimulation  
    from esys.geometry import RectangularDomain,VectorConstrainer  
    from esys.probe import Probe  
    from esys.input import InterpolatedTimeProfile,GausseanProfile  
   
    dom=RectangularDomain()  
    constraints=VectorConstrainer()  
    constraints.domain=Link(dom)  
    constraints.left=[1,1,1]  
    constraints.right=[1,1,1]  
    constraints.top=[1,1,1]  
    constraints.bottom=[1,1,1]  
    constraints.front=[1,1,1]  
    constraints.back=[1,1,1]  
   
    flow=SteadyIncompressibleFlow()  
    flow.domain=Link(dom)  
    flow.velocity=0.  
    flow.pressure=0.  
    flow.viscocity=1.  
    flow.internal_force=[1.,1.,0.]  
    flow.location_prescribed_velocities=Link(constraints,"location_of_constraint")  
    flow.prescribed_velocities=Data(), \  
   
    ptest=Probe()  
    ptest.reference("x[0]+x[1]")  
    ptest.value=Link(flow,"pressure")  
   
    s=ExplicitSimulation([dom,constraints,flow,ptest])  
    s.writeXML()  
    s.run()  

Legend:
Removed from v.144  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26