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

Contents of /trunk/esys2/escript/py_src/modellib/flow.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 142 - (show annotations)
Mon Jul 25 05:28:20 2005 UTC (14 years ago) by jgs
File MIME type: text/x-python
File size: 4785 byte(s)
Merge of development branch back to main trunk on 2005-07-25

1 # $Id$
2
3
4 from esys.escript import *
5 from esys.modelframe import Model,IterationDivergenceError
6 from esys.linearPDEs import LameEquation
7
8 class SteadyIncompressibleFlow(Model):
9 """
10
11 \f[
12 -\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i
13 \f]
14
15 \f[
16 \sigma_{ij}=2\eta D_{ij}-p\,\delta_{ij}
17 \f[
18 D_{ij}=\frac{1}{2}\left( v_{j,i} + v_{i,j }\right) \; .
19 \f]
20
21 \f[
22 -v_{k,k} = 0 \; .
23 \f]
24
25 """
26
27 def __init__(self,debug=False):
28 """set up model"""
29 Model.__init__(self,debug=debug)
30 self.declareParameter(domain=None, \
31 velocity=0., \
32 pressure=0., \
33 viscocity=1., \
34 internal_force=0., \
35 location_prescribed_velocities=Data(), \
36 prescribed_velocities=Data(), \
37 rel_tol=1.e-8,abs_tol=1.e-15,max_iter=10,relaxation=0.01)
38 self.__iter=0
39
40 def doInitialization(self,t):
41 """initialize model"""
42 self.__pressure_history=None
43 self.__dt_history=None
44 self.__pde=LameEquation(self.domain)
45
46 def stress(self):
47 """returns current stress"""
48 return 2*self.viscocity*self.stretching-self.pressure
49
50 def stretching(self):
51 """returns stertching tensor"""
52 g=grad(self.velocity)
53 return (g+transpose(g))/2
54
55 def doIterationInitialization(self,dt):
56 """
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 get an initial guess
61 """
62 if self.__dt_history==None:
63 self.__pressure_history=self.pressure
64 self.__iter=0
65 else:
66 self.__pressure_history,self.pressure=self.pressure,(1.+dt/self.__dt_history)*self.pressure+dt/self.__dt_history*self.__pressure_history
67 self.__iter=1
68 self.diff=1.e400
69
70 def doIterationStep(self,dt):
71 """
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 self.__pde.setValue(lame_lambda=self.viscosity, \
79 lame_my=penalty, \
80 F=F, \
81 sigma0=self.pressure*kronecker(self.__pde.getDomain()), \
82 r=self.location_prescribed_velocities, \
83 q=self.location_prescribed_velocities)
84 self.__pde.setTolerance(self.tol*1.e-2)
85 self.velocity=self.pde.getSolution()
86 update=penalty*div(self.velocity)
87 self.pressure=self.pressure-update
88 self.diff,diff=Lsup(update),self.diff
89 print "Pressure iteration: step %d: correction %e"%(self.__iter,self.diff/Lsup(self.pressure))
90 if diff<=self.diff:
91 raise IterationDivergenceError,"no improvement in pressure iteration"
92
93 def terminate(self):
94 """iteration is terminated if relative pressure change is less then rel_tol"""
95 if self.iter<2:
96 return False
97 else:
98 return self.diff<self.rel_tol*Lsup(self.pressure)+self.abs_tol
99
100 def doIterationFinalization(self,dt):
101 self.__dt_history=dt
102 if __name__=="__main__":
103
104
105
106 from esys.modelframe import Link,Simulation,ExplicitSimulation
107 from esys.geometry import RectangularDomain,VectorConstrainer
108 from esys.probe import Probe
109 from esys.input import InterpolatedTimeProfile,GausseanProfile
110
111 dom=RectangularDomain()
112 constraints=VectorConstrainer()
113 constraints.domain=Link(dom)
114 constraints.left=[1,1,1]
115 constraints.right=[1,1,1]
116 constraints.top=[1,1,1]
117 constraints.bottom=[1,1,1]
118 constraints.front=[1,1,1]
119 constraints.back=[1,1,1]
120
121 flow=SteadyIncompressibleFlow()
122 flow.domain=Link(dom)
123 flow.velocity=0.
124 flow.pressure=0.
125 flow.viscocity=1.
126 flow.internal_force=[1.,1.,0.]
127 flow.location_prescribed_velocities=Link(constraints,"location_of_constraint")
128 flow.prescribed_velocities=Data(), \
129
130 ptest=Probe()
131 ptest.reference("x[0]+x[1]")
132 ptest.value=Link(flow,"pressure")
133
134 s=ExplicitSimulation([dom,constraints,flow,ptest])
135 s.writeXML()
136 s.run()

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26