/[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 2625 - (hide annotations)
Fri Aug 21 06:30:25 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 5089 byte(s)
Merging changes from new doco branch.
All docstrings are now in reStructured text.
A few email addresses have been fixes as well.
1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 jgs 144
14 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18 elspeth 628 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 elspeth 628
22 jgs 149 from esys.escript import *
23     from esys.escript.modelframe import Model,IterationDivergenceError
24     from esys.escript.linearPDEs import LameEquation
25 jgs 144
26     class SteadyIncompressibleFlow(Model):
27 jgs 147 """
28 jgs 144
29 jfenwick 2625 *-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i*
30 jgs 149
31 jfenwick 2625 *\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}*
32 jgs 144
33 jfenwick 2625 *D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right)*
34 jgs 149
35 jfenwick 2625 *v_{k,k} = 0*
36 jgs 144
37     """
38    
39 gross 918 def __init__(self,**kwargs):
40 jgs 149 """
41     set up model
42     """
43 gross 918 Model.__init__(self,**kwargs)
44 jgs 144 self.declareParameter(domain=None, \
45     velocity=0., \
46     pressure=0., \
47 jgs 147 viscosity=1., \
48 jgs 144 internal_force=0., \
49 gross 918 location_prescribed_velocity=None, \
50     prescribed_velocity=None, \
51 gross 820 rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.0001)
52 jgs 144 self.__iter=0
53    
54 jgs 147 def doInitialization(self):
55 jgs 149 """
56     initialize model
57     """
58 jgs 147 self.__p_old=None
59     self.__p_very_old=None
60     self.__dt_old=None
61 jgs 144 self.__pde=LameEquation(self.domain)
62 gross 2474 self.__pde.getSolverOptions().setSolverMethod(self.__pde.getSolverOptions().DIRECT)
63 gross 918 if self.location_prescribed_velocity == None: self.location_prescribed_velocit=Data()
64     if self.prescribed_velocity == None: self.prescribed_velocity=Data()
65 jgs 144
66     def stress(self):
67 jgs 149 """
68     returns current stress"""
69 jgs 147 return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain)
70 jgs 144
71     def stretching(self):
72 jgs 149 """
73     returns stertching tensor
74     """
75 jgs 144 g=grad(self.velocity)
76     return (g+transpose(g))/2
77 jgs 147
78     def doStepPreprocessing(self,dt):
79 jgs 144 """
80 jgs 149 step up pressure iteration
81 jgs 144
82 jgs 149 if run within a time dependend problem extrapolation of pressure from previous time steps is used to
83     get an initial guess (that needs some work!!!!!!!)
84 jgs 144 """
85 jgs 147 self.__iter=0
86     self.__diff=1.e40
87     if not self.__p_old==None:
88     if self.__p_very_old==None:
89     self.pressure=self.__p_old
90     else:
91     self.pressure=(1.+dt/self.__dt_old)*self.__p_old-dt/self.__dt_old*self.__p_very_old
92 jgs 144
93 jgs 147 def doStep(self,dt):
94 jgs 144 """
95    
96 jgs 149 performs an iteration step of the penalty method.
97     IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached.
98 jgs 144
99     """
100     penalty=self.viscosity/self.relaxation
101 jgs 147 self.__pde.setValue(lame_mu=self.viscosity, \
102     lame_lambda=penalty, \
103     F=self.internal_force, \
104     sigma=self.pressure*kronecker(self.__pde.getDomain()), \
105     r=self.prescribed_velocity, \
106     q=self.location_prescribed_velocity)
107 gross 2474 self.__pde.getSolverOptions().setTolerance(self.rel_tol/10.)
108 jgs 147 self.velocity=self.__pde.getSolution()
109 jgs 144 update=penalty*div(self.velocity)
110     self.pressure=self.pressure-update
111 jgs 147 self.__diff,diff_old=Lsup(update),self.__diff
112     self.__iter+=1
113     self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
114     self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
115     self.trace("pressure correction: %e"%self.__diff)
116     if self.__iter>2 and diff_old<self.__diff:
117     self.trace("Pressure iteration failed!")
118 jgs 144 raise IterationDivergenceError,"no improvement in pressure iteration"
119 jgs 147 if self.__iter>self.max_iter:
120     raise IterationDivergenceError,"Maximum number of iterations steps reached"
121 jgs 144
122 jgs 147 def terminateIteration(self):
123     """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
124     return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol
125 jgs 144
126 jgs 147 def doStepPostprocessing(self,dt):
127     self.__dt_old=dt
128     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