/[escript]/branches/3.4.1/modellib/py_src/flow.py
ViewVC logotype

Annotation of /branches/3.4.1/modellib/py_src/flow.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4587 - (hide annotations)
Wed Dec 11 06:17:09 2013 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 5188 byte(s)
Preparation begins

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