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	# $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