/[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 144 - (hide annotations)
Mon Jul 25 05:42:21 2005 UTC (15 years, 4 months ago) by jgs
Original Path: trunk/esys2/modellib/py_src/flow.py
File MIME type: text/x-python
File size: 4785 byte(s)
moved form old modellib directory

1 jgs 144 # $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