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

Annotation of /trunk/modellib/py_src/flow.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (hide annotations)
Tue Aug 23 01:24:31 2005 UTC (15 years, 3 months ago) by jgs
Original Path: trunk/esys2/modellib/py_src/flow.py
File MIME type: text/x-python
File size: 4176 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

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 2*self.viscosity*self.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_old*self.__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     self.pressure=self.pressure-update
88 jgs 147 self.__diff,diff_old=Lsup(update),self.__diff
89     self.__iter+=1
90     self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
91     self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
92     self.trace("pressure correction: %e"%self.__diff)
93     if self.__iter>2 and diff_old<self.__diff:
94     self.trace("Pressure iteration failed!")
95 jgs 144 raise IterationDivergenceError,"no improvement in pressure iteration"
96 jgs 147 if self.__iter>self.max_iter:
97     raise IterationDivergenceError,"Maximum number of iterations steps reached"
98 jgs 144
99 jgs 147 def terminateIteration(self):
100     """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
101     return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol
102 jgs 144
103 jgs 147 def doStepPostprocessing(self,dt):
104     self.__dt_old=dt
105     self.__p_old,self.__p_very_old=self.pressure,self.__p_old

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26