modellib/py_src/flow.py

# $Id$

from escript.escript import *
from escript.modelframe import Model,IterationDivergenceError
from escript.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., \
                                 viscosity=1., \
                                 internal_force=0., \
                                 location_prescribed_velocity=Data(), \
                                 prescribed_velocity=Data(), \
                                 rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.001)
           self.__iter=0

       def doInitialization(self):
           """initialize model"""
           self.__p_old=None
           self.__p_very_old=None
           self.__dt_old=None
           self.__pde=LameEquation(self.domain)

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

       def stretching(self):
           """returns stertching tensor"""
           g=grad(self.velocity)
           return (g+transpose(g))/2

       def doStepPreprocessing(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 (that needs some work!!!!!!!)
            """
            self.__iter=0
            self.__diff=1.e40
            if not self.__p_old==None:
               if self.__p_very_old==None:
                  self.pressure=self.__p_old
               else:
                  self.pressure=(1.+dt/self.__dt_old)*self.__p_old-dt/self.__dt_old*self.__p_very_old

       def doStep(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_mu=self.viscosity, \
                              lame_lambda=penalty, \
                              F=self.internal_force, \
                              sigma=self.pressure*kronecker(self.__pde.getDomain()), \
                              r=self.prescribed_velocity, \
                              q=self.location_prescribed_velocity)
          self.__pde.setTolerance(self.rel_tol/10.)
          self.velocity=self.__pde.getSolution()
          update=penalty*div(self.velocity)
          print "f=",inf(update),sup(update)
          self.pressure=self.pressure-update
          self.__diff,diff_old=Lsup(update),self.__diff
          self.__iter+=1
          self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
          self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
          self.trace("pressure correction: %e"%self.__diff)
          if self.__iter>2 and diff_old<self.__diff:
              self.trace("Pressure iteration failed!")
              raise IterationDivergenceError,"no improvement in pressure iteration"
          if self.__iter>self.max_iter:
              raise IterationDivergenceError,"Maximum number of iterations steps reached"

       def terminateIteration(self):
          """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
          return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol

       def doStepPostprocessing(self,dt):
          self.__dt_old=dt
          self.__p_old,self.__p_very_old=self.pressure,self.__p_old
1	jgs	144	# $Id$
2
3	jgs	147	from escript.escript import *
4			from escript.modelframe import Model,IterationDivergenceError
5			from escript.linearPDEs import LameEquation
6	jgs	144
7			class SteadyIncompressibleFlow(Model):
8	jgs	147	"""
9	jgs	144
10			\f[
11			-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i
12			\f]
13
14			\f[
15			\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}
16			\f[
17			D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; .
18			\f]
19
20			\f[
21			-v_{k,k} = 0 \; .
22			\f]
23	jgs	147
24	jgs	144	"""
25
26			def __init__(self,debug=False):
27			"""set up model"""
28			Model.__init__(self,debug=debug)
29			self.declareParameter(domain=None, \
30			velocity=0., \
31			pressure=0., \
32	jgs	147	viscosity=1., \
33	jgs	144	internal_force=0., \
34	jgs	147	location_prescribed_velocity=Data(), \
35			prescribed_velocity=Data(), \
36			rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.001)
37	jgs	144	self.__iter=0
38
39	jgs	147	def doInitialization(self):
40	jgs	144	"""initialize model"""
41	jgs	147	self.__p_old=None
42			self.__p_very_old=None
43			self.__dt_old=None
44	jgs	144	self.__pde=LameEquation(self.domain)
45
46			def stress(self):
47			"""returns current stress"""
48	jgs	147	return 2self.viscosityself.stretching-self.pressure*kronecker(self.domain)
49	jgs	144
50			def stretching(self):
51			"""returns stertching tensor"""
52			g=grad(self.velocity)
53			return (g+transpose(g))/2
54	jgs	147
55			def doStepPreprocessing(self,dt):
56	jgs	144	"""
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	jgs	147	get an initial guess (that needs some work!!!!!!!)
61	jgs	144	"""
62	jgs	147	self.__iter=0
63			self.__diff=1.e40
64			if not self.__p_old==None:
65			if self.__p_very_old==None:
66			self.pressure=self.__p_old
67			else:
68			self.pressure=(1.+dt/self.__dt_old)self.__p_old-dt/self.__dt_oldself.__p_very_old
69	jgs	144
70	jgs	147	def doStep(self,dt):
71	jgs	144	"""
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	jgs	147	self.__pde.setValue(lame_mu=self.viscosity, \
79			lame_lambda=penalty, \
80			F=self.internal_force, \
81			sigma=self.pressure*kronecker(self.__pde.getDomain()), \
82			r=self.prescribed_velocity, \
83			q=self.location_prescribed_velocity)
84			self.__pde.setTolerance(self.rel_tol/10.)
85			self.velocity=self.__pde.getSolution()
86	jgs	144	update=penalty*div(self.velocity)
87	jgs	147	print "f=",inf(update),sup(update)
88	jgs	144	self.pressure=self.pressure-update
89	jgs	147	self.__diff,diff_old=Lsup(update),self.__diff
90			self.__iter+=1
91			self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
92			self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
93			self.trace("pressure correction: %e"%self.__diff)
94			if self.__iter>2 and diff_old<self.__diff:
95			self.trace("Pressure iteration failed!")
96	jgs	144	raise IterationDivergenceError,"no improvement in pressure iteration"
97	jgs	147	if self.__iter>self.max_iter:
98			raise IterationDivergenceError,"Maximum number of iterations steps reached"
99	jgs	144
100	jgs	147	def terminateIteration(self):
101			"""iteration is terminateIterationd if relative pressure change is less then rel_tol"""
102			return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol
103	jgs	144
104	jgs	147	def doStepPostprocessing(self,dt):
105			self.__dt_old=dt
106			self.__p_old,self.__p_very_old=self.pressure,self.__p_old
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision