/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

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

revision 1481 by artak, Wed Apr 9 00:45:47 2008 UTC revision 1841 by gross, Fri Oct 3 03:57:52 2008 UTC
# Line 1  Line 1 
1  # $Id:$  
2  #  ########################################################
 #######################################################  
 #  
 #       Copyright 2008 by University of Queensland  
3  #  #
4  #                http://esscc.uq.edu.au  # Copyright (c) 2003-2008 by University of Queensland
5  #        Primary Business: Queensland, Australia  # Earth Systems Science Computational Center (ESSCC)
6  #  Licensed under the Open Software License version 3.0  # http://www.uq.edu.au/esscc
 #     http://www.opensource.org/licenses/osl-3.0.php  
7  #  #
8  #######################################################  # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11  #  #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 31  Some models for flow
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
34    
35  from escript import *  from escript import *
36  import util  import util
37  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
38  from pdetools import HomogeneousSaddlePointProblem  from pdetools import HomogeneousSaddlePointProblem,Projector
39    
40    class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):
41          """
42          solves
43    
44              -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
45                    u_{i,i}=0
46    
47              u=0 where  fixed_u_mask>0
48              eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
49    
50          if surface_stress is not give 0 is assumed.
51    
52          typical usage:
53    
54                sp=StokesProblemCartesian(domain)
55                sp.setTolerance()
56                sp.initialize(...)
57                v,p=sp.solve(v0,p0)
58          """
59          def __init__(self,domain,**kwargs):
60             HomogeneousSaddlePointProblem.__init__(self,**kwargs)
61             self.domain=domain
62             self.vol=util.integrate(1.,Function(self.domain))
63             self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
64             self.__pde_u.setSymmetryOn()
65             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
66                
67             # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
68             # self.__pde_proj.setReducedOrderOn()
69             # self.__pde_proj.setSymmetryOn()
70             # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)
71    
72          def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
73            self.eta=eta
74            A =self.__pde_u.createCoefficient("A")
75        self.__pde_u.setValue(A=Data())
76            for i in range(self.domain.getDim()):
77            for j in range(self.domain.getDim()):
78                A[i,j,j,i] += 1.
79                A[i,j,i,j] += 1.
80            # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
81            self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
82    
83            # self.__pde_proj.setValue(D=1/eta)
84            # self.__pde_proj.setValue(Y=1.)
85            # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))
86            self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
87    
88          def B(self,arg):
89             a=util.div(arg, ReducedFunction(self.domain))
90             return a-util.integrate(a)/self.vol
91    
92          def inner(self,p0,p1):
93             return util.integrate(p0*p1)
94             # return util.integrate(1/self.__inv_eta**2*p0*p1)
95    
96          def getStress(self,u):
97             mg=util.grad(u)
98             return 2.*self.eta*util.symmetric(mg)
99          def getEtaEffective(self):
100             return self.eta
101    
102          def solve_A(self,u,p):
103             """
104             solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
105             """
106             self.__pde_u.setTolerance(self.getSubProblemTolerance())
107             self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))
108             return  self.__pde_u.getSolution(verbose=self.show_details)
109    
110    
111          def solve_prec(self,p):
112            a=self.__inv_eta*p
113            return a-util.integrate(a)/self.vol
114    
115          def stoppingcriterium(self,Bv,v,p):
116              n_r=util.sqrt(self.inner(Bv,Bv))
117              n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
118              if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)
119              if self.iter == 0: self.__n_v=n_v;
120              self.__n_v, n_v_old =n_v, self.__n_v
121              self.iter+=1
122              if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
123                  if self.verbose: print "PCG terminated after %s steps."%self.iter
124                  return True
125              else:
126                  return False
127          def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
128          if TOL==None:
129                 TOL=self.getTolerance()
130              if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
131              self.iter+=1
132              
133              if norm_r <= norm_b*TOL:
134                  if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
135                  return True
136              else:
137                  return False
138    
139    
140  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
141        """        """
# Line 63  class StokesProblemCartesian(Homogeneous Line 162  class StokesProblemCartesian(Homogeneous
162           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
163           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
164           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
165           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
166                            
167           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
168           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 76  class StokesProblemCartesian(Homogeneous Line 175  class StokesProblemCartesian(Homogeneous
175    
176        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
177          self.eta=eta          self.eta=eta
178          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
179      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
180          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
181          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
182              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
183              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
184      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
185          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
186    
187        def B(self,arg):        def B(self,arg):
# Line 96  class StokesProblemCartesian(Homogeneous Line 195  class StokesProblemCartesian(Homogeneous
195           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
196           return util.integrate(s0*s1)           return util.integrate(s0*s1)
197    
198          def inner_a(self,a0,a1):
199             p0=util.interpolate(a0[1],Function(self.domain))
200             p1=util.interpolate(a1[1],Function(self.domain))
201             alfa=(1/self.vol)*util.integrate(p0)
202             beta=(1/self.vol)*util.integrate(p1)
203         v0=util.grad(a0[0])
204         v1=util.grad(a1[0])
205             return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
206    
207    
208        def getStress(self,u):        def getStress(self,u):
209           mg=util.grad(u)           mg=util.grad(u)
210           return 2.*self.eta*util.symmetric(mg)           return 2.*self.eta*util.symmetric(mg)
211          def getEtaEffective(self):
212             return self.eta
213    
214        def solve_A(self,u,p):        def solve_A(self,u,p):
215           """           """
# Line 108  class StokesProblemCartesian(Homogeneous Line 219  class StokesProblemCartesian(Homogeneous
219           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
220           return  self.__pde_u.getSolution(verbose=self.show_details)           return  self.__pde_u.getSolution(verbose=self.show_details)
221    
222    
223        def solve_prec(self,p):        def solve_prec(self,p):
224         #proj=Projector(domain=self.domain, reduce = True, fast=False)
225           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
226           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
227           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
228           return q           return q
229    
230          def solve_prec1(self,p):
231         #proj=Projector(domain=self.domain, reduce = True, fast=False)
232             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
233             self.__pde_prec.setValue(Y=p)
234             q=self.__pde_prec.getSolution(verbose=self.show_details)
235         q0=util.interpolate(q,Function(self.domain))
236             print util.inf(q*q0),util.sup(q*q0)
237             q-=(1/self.vol)*util.integrate(q0)
238             print util.inf(q*q0),util.sup(q*q0)
239             return q
240    
241        def stoppingcriterium(self,Bv,v,p):        def stoppingcriterium(self,Bv,v,p):
242            n_r=util.sqrt(self.inner(Bv,Bv))            n_r=util.sqrt(self.inner(Bv,Bv))
243            n_v=util.Lsup(v)            n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
244            if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)            if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
245              if self.iter == 0: self.__n_v=n_v;
246              self.__n_v, n_v_old =n_v, self.__n_v
247            self.iter+=1            self.iter+=1
248            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
249                if self.verbose: print "PCG terminated after %s steps."%self.iter                if self.verbose: print "PCG terminated after %s steps."%self.iter
250                return True                return True
251            else:            else:
252                return False                return False
253        def stoppingcriterium_GMRES(self,norm_r,norm_b):        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
254            if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())        if TOL==None:
255                 TOL=self.getTolerance()
256              if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
257            self.iter+=1            self.iter+=1
258            if norm_r <= norm_b*self.getTolerance():            
259                if self.verbose: print "GMRES terminated after %s steps."%self.iter            if norm_r <= norm_b*TOL:
260                  if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
261                return True                return True
262            else:            else:
263                return False                return False
264    
265        def stoppingcriterium_MINRES(self,norm_r,norm_Ax):  
           if self.verbose: print "MINRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_Ax*self.getTolerance())  
           self.iter+=1  
           if norm_r <= norm_Ax*self.getTolerance():  
               if self.verbose: print "MINRES terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  

Legend:
Removed from v.1481  
changed lines
  Added in v.1841

  ViewVC Help
Powered by ViewVC 1.1.26