/[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 147 - (hide annotations)
Fri Aug 12 01:45:47 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: 4221 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26