modellib/py_src/flow.py

# $Id$


from esys.escript import *
from esys.modelframe import Model,IterationDivergenceError
from esys.linearPDEs import LameEquation

class SteadyIncompressibleFlow(Model):
       """ 

                   \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]
                   
       """

       def __init__(self,debug=False):
           """set up model"""
           Model.__init__(self,debug=debug)
           self.declareParameter(domain=None, \
                                 velocity=0., \
                                 pressure=0., \
                                 viscocity=1., \
                                 internal_force=0., \
                                 location_prescribed_velocities=Data(), \
                                 prescribed_velocities=Data(), \
                                 rel_tol=1.e-8,abs_tol=1.e-15,max_iter=10,relaxation=0.01) 
           self.__iter=0

       def doInitialization(self,t):
           """initialize model"""
           self.__pressure_history=None
           self.__dt_history=None
           self.__pde=LameEquation(self.domain)

       def stress(self):
           """returns current stress"""
           return 2*self.viscocity*self.stretching-self.pressure

       def stretching(self):
           """returns stertching tensor"""
           g=grad(self.velocity)
           return (g+transpose(g))/2
           
       def doIterationInitialization(self,dt):
            """
               step up pressure iteration

               if run within a time dependend problem extrapolation of pressure from previous time steps is used to
               get an initial guess 
            """
            if self.__dt_history==None:
               self.__pressure_history=self.pressure
               self.__iter=0
            else:
               self.__pressure_history,self.pressure=self.pressure,(1.+dt/self.__dt_history)*self.pressure+dt/self.__dt_history*self.__pressure_history
               self.__iter=1
            self.diff=1.e400

       def doIterationStep(self,dt):
          """

             performs an iteration step of the penalty method.
             IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.

          """
          penalty=self.viscosity/self.relaxation
          self.__pde.setValue(lame_lambda=self.viscosity, \
                              lame_my=penalty, \
                              F=F, \
                              sigma0=self.pressure*kronecker(self.__pde.getDomain()), \
                              r=self.location_prescribed_velocities, \
                              q=self.location_prescribed_velocities) 
          self.__pde.setTolerance(self.tol*1.e-2)
          self.velocity=self.pde.getSolution()
          update=penalty*div(self.velocity)
          self.pressure=self.pressure-update
          self.diff,diff=Lsup(update),self.diff
          print "Pressure iteration: step %d: correction %e"%(self.__iter,self.diff/Lsup(self.pressure))
          if diff<=self.diff:
              raise IterationDivergenceError,"no improvement in pressure iteration"

       def terminate(self):
          """iteration is terminated if relative pressure change is less then rel_tol"""
          if self.iter<2:
              return False
          else:
             return self.diff<self.rel_tol*Lsup(self.pressure)+self.abs_tol

       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()
1	jgs	144	# $Id$
2
3
4			from esys.escript import *
5			from esys.modelframe import Model,IterationDivergenceError
6			from esys.linearPDEs import LameEquation
7
8			class SteadyIncompressibleFlow(Model):
9			"""
10
11			\f[
12			-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i
13			\f]
14
15			\f[
16			\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}
17			\f[
18			D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; .
19			\f]
20
21			\f[
22			-v_{k,k} = 0 \; .
23			\f]
24
25			"""
26
27			def __init__(self,debug=False):
28			"""set up model"""
29			Model.__init__(self,debug=debug)
30			self.declareParameter(domain=None, \
31			velocity=0., \
32			pressure=0., \
33			viscocity=1., \
34			internal_force=0., \
35			location_prescribed_velocities=Data(), \
36			prescribed_velocities=Data(), \
37			rel_tol=1.e-8,abs_tol=1.e-15,max_iter=10,relaxation=0.01)
38			self.__iter=0
39
40			def doInitialization(self,t):
41			"""initialize model"""
42			self.__pressure_history=None
43			self.__dt_history=None
44			self.__pde=LameEquation(self.domain)
45
46			def stress(self):
47			"""returns current stress"""
48			return 2self.viscocityself.stretching-self.pressure
49
50			def stretching(self):
51			"""returns stertching tensor"""
52			g=grad(self.velocity)
53			return (g+transpose(g))/2
54
55			def doIterationInitialization(self,dt):
56			"""
57			step up pressure iteration
58
59			if run within a time dependend problem extrapolation of pressure from previous time steps is used to
60			get an initial guess
61			"""
62			if self.__dt_history==None:
63			self.__pressure_history=self.pressure
64			self.__iter=0
65			else:
66			self.__pressure_history,self.pressure=self.pressure,(1.+dt/self.__dt_history)self.pressure+dt/self.__dt_historyself.__pressure_history
67			self.__iter=1
68			self.diff=1.e400
69
70			def doIterationStep(self,dt):
71			"""
72
73			performs an iteration step of the penalty method.
74			IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
75
76			"""
77			penalty=self.viscosity/self.relaxation
78			self.__pde.setValue(lame_lambda=self.viscosity, \
79			lame_my=penalty, \
80			F=F, \
81			sigma0=self.pressure*kronecker(self.__pde.getDomain()), \
82			r=self.location_prescribed_velocities, \
83			q=self.location_prescribed_velocities)
84			self.__pde.setTolerance(self.tol*1.e-2)
85			self.velocity=self.pde.getSolution()
86			update=penalty*div(self.velocity)
87			self.pressure=self.pressure-update
88			self.diff,diff=Lsup(update),self.diff
89			print "Pressure iteration: step %d: correction %e"%(self.__iter,self.diff/Lsup(self.pressure))
90			if diff<=self.diff:
91			raise IterationDivergenceError,"no improvement in pressure iteration"
92
93			def terminate(self):
94			"""iteration is terminated if relative pressure change is less then rel_tol"""
95			if self.iter<2:
96			return False
97			else:
98			return self.diff<self.rel_tol*Lsup(self.pressure)+self.abs_tol
99
100			def doIterationFinalization(self,dt):
101			self.__dt_history=dt
102			if __name__=="__main__":
103
104
105
106			from esys.modelframe import Link,Simulation,ExplicitSimulation
107			from esys.geometry import RectangularDomain,VectorConstrainer
108			from esys.probe import Probe
109			from esys.input import InterpolatedTimeProfile,GausseanProfile
110
111			dom=RectangularDomain()
112			constraints=VectorConstrainer()
113			constraints.domain=Link(dom)
114			constraints.left=[1,1,1]
115			constraints.right=[1,1,1]
116			constraints.top=[1,1,1]
117			constraints.bottom=[1,1,1]
118			constraints.front=[1,1,1]
119			constraints.back=[1,1,1]
120
121			flow=SteadyIncompressibleFlow()
122			flow.domain=Link(dom)
123			flow.velocity=0.
124			flow.pressure=0.
125			flow.viscocity=1.
126			flow.internal_force=[1.,1.,0.]
127			flow.location_prescribed_velocities=Link(constraints,"location_of_constraint")
128			flow.prescribed_velocities=Data(), \
129
130			ptest=Probe()
131			ptest.reference("x[0]+x[1]")
132			ptest.value=Link(flow,"pressure")
133
134			s=ExplicitSimulation([dom,constraints,flow,ptest])
135			s.writeXML()
136			s.run()
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision