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

Diff of /trunk/modellib/py_src/flow.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 144 by jgs, Mon Jul 25 05:42:21 2005 UTC revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
2    
3    from escript.escript import *
4  from esys.escript import *  from escript.modelframe import Model,IterationDivergenceError
5  from esys.modelframe import Model,IterationDivergenceError  from escript.linearPDEs import LameEquation
 from esys.linearPDEs import LameEquation  
6    
7  class SteadyIncompressibleFlow(Model):  class SteadyIncompressibleFlow(Model):
8         """         """
9    
10                     \f[                     \f[
11                     -\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i                     -\left(\eta\left(v_{i,j}+v_{j,i}\right)\right)_{,j}+p_{,i}=F_i
# Line 21  class SteadyIncompressibleFlow(Model): Line 20  class SteadyIncompressibleFlow(Model):
20                     \f[                     \f[
21                     -v_{k,k} = 0 \; .                     -v_{k,k} = 0 \; .
22                     \f]                     \f]
23                      
24         """         """
25    
26         def __init__(self,debug=False):         def __init__(self,debug=False):
# Line 30  class SteadyIncompressibleFlow(Model): Line 29  class SteadyIncompressibleFlow(Model):
29             self.declareParameter(domain=None, \             self.declareParameter(domain=None, \
30                                   velocity=0., \                                   velocity=0., \
31                                   pressure=0., \                                   pressure=0., \
32                                   viscocity=1., \                                   viscosity=1., \
33                                   internal_force=0., \                                   internal_force=0., \
34                                   location_prescribed_velocities=Data(), \                                   location_prescribed_velocity=Data(), \
35                                   prescribed_velocities=Data(), \                                   prescribed_velocity=Data(), \
36                                   rel_tol=1.e-8,abs_tol=1.e-15,max_iter=10,relaxation=0.01)                                   rel_tol=1.e-3,abs_tol=0.,max_iter=10,relaxation=0.001)
37             self.__iter=0             self.__iter=0
38    
39         def doInitialization(self,t):         def doInitialization(self):
40             """initialize model"""             """initialize model"""
41             self.__pressure_history=None             self.__p_old=None
42             self.__dt_history=None             self.__p_very_old=None
43               self.__dt_old=None
44             self.__pde=LameEquation(self.domain)             self.__pde=LameEquation(self.domain)
45    
46         def stress(self):         def stress(self):
47             """returns current stress"""             """returns current stress"""
48             return 2*self.viscocity*self.stretching-self.pressure             return 2*self.viscosity*self.stretching-self.pressure*kronecker(self.domain)
49    
50         def stretching(self):         def stretching(self):
51             """returns stertching tensor"""             """returns stertching tensor"""
52             g=grad(self.velocity)             g=grad(self.velocity)
53             return (g+transpose(g))/2             return (g+transpose(g))/2
54              
55         def doIterationInitialization(self,dt):         def doStepPreprocessing(self,dt):
56              """              """
57                 step up pressure iteration                 step up pressure iteration
58    
59                 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
60                 get an initial guess                 get an initial guess (that needs some work!!!!!!!)
61              """              """
62              if self.__dt_history==None:              self.__iter=0
63                 self.__pressure_history=self.pressure              self.__diff=1.e40
64                 self.__iter=0              if not self.__p_old==None:
65              else:                 if self.__p_very_old==None:
66                 self.__pressure_history,self.pressure=self.pressure,(1.+dt/self.__dt_history)*self.pressure+dt/self.__dt_history*self.__pressure_history                    self.pressure=self.__p_old
67                 self.__iter=1                 else:
68              self.diff=1.e400                    self.pressure=(1.+dt/self.__dt_old)*self.__p_old-dt/self.__dt_old*self.__p_very_old
69    
70         def doIterationStep(self,dt):         def doStep(self,dt):
71            """            """
72    
73               performs an iteration step of the penalty method.               performs an iteration step of the penalty method.
# Line 75  class SteadyIncompressibleFlow(Model): Line 75  class SteadyIncompressibleFlow(Model):
75    
76            """            """
77            penalty=self.viscosity/self.relaxation            penalty=self.viscosity/self.relaxation
78            self.__pde.setValue(lame_lambda=self.viscosity, \            self.__pde.setValue(lame_mu=self.viscosity, \
79                                lame_my=penalty, \                                lame_lambda=penalty, \
80                                F=F, \                                F=self.internal_force, \
81                                sigma0=self.pressure*kronecker(self.__pde.getDomain()), \                                sigma=self.pressure*kronecker(self.__pde.getDomain()), \
82                                r=self.location_prescribed_velocities, \                                r=self.prescribed_velocity, \
83                                q=self.location_prescribed_velocities)                                q=self.location_prescribed_velocity)
84            self.__pde.setTolerance(self.tol*1.e-2)            self.__pde.setTolerance(self.rel_tol/10.)
85            self.velocity=self.pde.getSolution()            self.velocity=self.__pde.getSolution()
86            update=penalty*div(self.velocity)            update=penalty*div(self.velocity)
87              print "f=",inf(update),sup(update)
88            self.pressure=self.pressure-update            self.pressure=self.pressure-update
89            self.diff,diff=Lsup(update),self.diff            self.__diff,diff_old=Lsup(update),self.__diff
90            print "Pressure iteration: step %d: correction %e"%(self.__iter,self.diff/Lsup(self.pressure))            self.__iter+=1
91            if diff<=self.diff:            self.trace("velocity range %e:%e"%(inf(self.velocity),sup(self.velocity)))
92              self.trace("pressure range %e:%e"%(inf(self.pressure),sup(self.pressure)))
93              self.trace("pressure correction: %e"%self.__diff)
94              if self.__iter>2 and diff_old<self.__diff:
95                  self.trace("Pressure iteration failed!")
96                raise IterationDivergenceError,"no improvement in pressure iteration"                raise IterationDivergenceError,"no improvement in pressure iteration"
97              if self.__iter>self.max_iter:
98                  raise IterationDivergenceError,"Maximum number of iterations steps reached"
99    
100         def terminate(self):         def terminateIteration(self):
101            """iteration is terminated if relative pressure change is less then rel_tol"""            """iteration is terminateIterationd if relative pressure change is less then rel_tol"""
102            if self.iter<2:            return self.__diff<=self.rel_tol*Lsup(self.pressure)+self.abs_tol
103                return False  
104            else:         def doStepPostprocessing(self,dt):
105               return self.diff<self.rel_tol*Lsup(self.pressure)+self.abs_tol            self.__dt_old=dt
106              self.__p_old,self.__p_very_old=self.pressure,self.__p_old
        def doIterationFinalization(self,dt):  
           self.__dt_history=dt  
 if __name__=="__main__":  
   
   
   
    from esys.modelframe import Link,Simulation,ExplicitSimulation  
    from esys.geometry import RectangularDomain,VectorConstrainer  
    from esys.probe import Probe  
    from esys.input import InterpolatedTimeProfile,GausseanProfile  
   
    dom=RectangularDomain()  
    constraints=VectorConstrainer()  
    constraints.domain=Link(dom)  
    constraints.left=[1,1,1]  
    constraints.right=[1,1,1]  
    constraints.top=[1,1,1]  
    constraints.bottom=[1,1,1]  
    constraints.front=[1,1,1]  
    constraints.back=[1,1,1]  
   
    flow=SteadyIncompressibleFlow()  
    flow.domain=Link(dom)  
    flow.velocity=0.  
    flow.pressure=0.  
    flow.viscocity=1.  
    flow.internal_force=[1.,1.,0.]  
    flow.location_prescribed_velocities=Link(constraints,"location_of_constraint")  
    flow.prescribed_velocities=Data(), \  
   
    ptest=Probe()  
    ptest.reference("x[0]+x[1]")  
    ptest.value=Link(flow,"pressure")  
   
    s=ExplicitSimulation([dom,constraints,flow,ptest])  
    s.writeXML()  
    s.run()  
Legend:
Removed from v.144  
changed lines
  Added in v.147
  ViewVC Help
Powered by ViewVC 1.1.26