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

Revision 142 - (show annotations)
Mon Jul 25 05:28:20 2005 UTC (16 years, 11 months 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

Properties

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