1 |
# $Id$ |
# $Id$ |
2 |
|
|
3 |
from escript.escript import * |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
4 |
from escript.modelframe import Model,IterationDivergenceError |
http://www.access.edu.au |
5 |
from escript.linearPDEs import LameEquation |
Primary Business: Queensland, Australia""" |
6 |
|
__license__="""Licensed under the Open Software License version 3.0 |
7 |
|
http://www.opensource.org/licenses/osl-3.0.php""" |
8 |
|
|
9 |
|
from esys.escript import * |
10 |
|
from esys.escript.modelframe import Model,IterationDivergenceError |
11 |
|
from esys.escript.linearPDEs import LameEquation |
12 |
|
|
13 |
class SteadyIncompressibleFlow(Model): |
class SteadyIncompressibleFlow(Model): |
14 |
""" |
""" |
15 |
|
|
16 |
\f[ |
M{-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i} |
17 |
-\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i |
|
18 |
\f] |
M{\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}} |
19 |
|
|
20 |
\f[ |
M{D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right)} |
21 |
\sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij} |
|
22 |
\f[ |
M{v_{k,k} = 0} |
|
D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; . |
|
|
\f] |
|
|
|
|
|
\f[ |
|
|
-v_{k,k} = 0 \; . |
|
|
\f] |
|
23 |
|
|
24 |
""" |
""" |
25 |
|
|
26 |
def __init__(self,debug=False): |
def __init__(self,**kwargs): |
27 |
"""set up model""" |
""" |
28 |
Model.__init__(self,debug=debug) |
set up model |
29 |
|
""" |
30 |
|
Model.__init__(self,**kwargs) |
31 |
self.declareParameter(domain=None, \ |
self.declareParameter(domain=None, \ |
32 |
velocity=0., \ |
velocity=0., \ |
33 |
pressure=0., \ |
pressure=0., \ |
34 |
viscosity=1., \ |
viscosity=1., \ |
35 |
internal_force=0., \ |
internal_force=0., \ |
36 |
location_prescribed_velocity=Data(), \ |
location_prescribed_velocity=None, \ |
37 |
prescribed_velocity=Data(), \ |
prescribed_velocity=None, \ |
38 |
rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.001) |
rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.0001) |
39 |
self.__iter=0 |
self.__iter=0 |
40 |
|
|
41 |
def doInitialization(self): |
def doInitialization(self): |
42 |
"""initialize model""" |
""" |
43 |
|
initialize model |
44 |
|
""" |
45 |
self.__p_old=None |
self.__p_old=None |
46 |
self.__p_very_old=None |
self.__p_very_old=None |
47 |
self.__dt_old=None |
self.__dt_old=None |
48 |
self.__pde=LameEquation(self.domain) |
self.__pde=LameEquation(self.domain) |
49 |
|
self.__pde.setSolverMethod(LameEquation.DIRECT) |
50 |
|
if self.location_prescribed_velocity == None: self.location_prescribed_velocit=Data() |
51 |
|
if self.prescribed_velocity == None: self.prescribed_velocity=Data() |
52 |
|
|
53 |
def stress(self): |
def stress(self): |
54 |
"""returns current stress""" |
""" |
55 |
|
returns current stress""" |
56 |
return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain) |
return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain) |
57 |
|
|
58 |
def stretching(self): |
def stretching(self): |
59 |
"""returns stertching tensor""" |
""" |
60 |
|
returns stertching tensor |
61 |
|
""" |
62 |
g=grad(self.velocity) |
g=grad(self.velocity) |
63 |
return (g+transpose(g))/2 |
return (g+transpose(g))/2 |
64 |
|
|
65 |
def doStepPreprocessing(self,dt): |
def doStepPreprocessing(self,dt): |
66 |
""" |
""" |
67 |
step up pressure iteration |
step up pressure iteration |
68 |
|
|
69 |
if run within a time dependend problem extrapolation of pressure from previous time steps is used to |
if run within a time dependend problem extrapolation of pressure from previous time steps is used to |
70 |
get an initial guess (that needs some work!!!!!!!) |
get an initial guess (that needs some work!!!!!!!) |
71 |
""" |
""" |
72 |
self.__iter=0 |
self.__iter=0 |
73 |
self.__diff=1.e40 |
self.__diff=1.e40 |
80 |
def doStep(self,dt): |
def doStep(self,dt): |
81 |
""" |
""" |
82 |
|
|
83 |
performs an iteration step of the penalty method. |
performs an iteration step of the penalty method. |
84 |
IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached. |
IterationDivergenceError is raised if pressure error cannot be reduced or max_iter is reached. |
85 |
|
|
86 |
""" |
""" |
87 |
penalty=self.viscosity/self.relaxation |
penalty=self.viscosity/self.relaxation |
94 |
self.__pde.setTolerance(self.rel_tol/10.) |
self.__pde.setTolerance(self.rel_tol/10.) |
95 |
self.velocity=self.__pde.getSolution() |
self.velocity=self.__pde.getSolution() |
96 |
update=penalty*div(self.velocity) |
update=penalty*div(self.velocity) |
|
print "f=",inf(update),sup(update) |
|
97 |
self.pressure=self.pressure-update |
self.pressure=self.pressure-update |
98 |
self.__diff,diff_old=Lsup(update),self.__diff |
self.__diff,diff_old=Lsup(update),self.__diff |
99 |
self.__iter+=1 |
self.__iter+=1 |