/[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 1673 by gross, Thu Jul 24 22:28:50 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC
# Line 1  Line 1 
1  # $Id:$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2008 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8  #  # Licensed under the Open Software License version 3.0
9  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
10  #  #
11    ########################################################
12    
13    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
16    Primary Business: Queensland, Australia"""
17    __license__="""Licensed under the Open Software License version 3.0
18    http://www.opensource.org/licenses/osl-3.0.php"""
19    __url__="http://www.uq.edu.au/esscc/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31    
32  __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:$"  
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
   
 class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  
       """  
       solves  
   
           -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i  
                 u_{i,i}=0  
   
           u=0 where  fixed_u_mask>0  
           eta*(u_{i,j}+u_{j,i})*n_j=surface_stress  
   
       if surface_stress is not give 0 is assumed.  
   
       typical usage:  
   
             sp=StokesProblemCartesian(domain)  
             sp.setTolerance()  
             sp.initialize(...)  
             v,p=sp.solve(v0,p0)  
       """  
       def __init__(self,domain,**kwargs):  
          HomogeneousSaddlePointProblem.__init__(self,**kwargs)  
          self.domain=domain  
          self.vol=util.integrate(1.,Function(self.domain))  
          self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())  
          self.__pde_u.setSymmetryOn()  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
          # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)  
          # self.__pde_proj.setReducedOrderOn()  
          # self.__pde_proj.setSymmetryOn()  
          # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)  
   
       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):  
         self.eta=eta  
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
         # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
   
         # self.__pde_proj.setValue(D=1/eta)  
         # self.__pde_proj.setValue(Y=1.)  
         # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))  
         self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
   
       def B(self,arg):  
          a=util.div(arg, ReducedFunction(self.domain))  
          return a-util.integrate(a)/self.vol  
   
       def inner(self,p0,p1):  
          return util.integrate(p0*p1)  
          # return util.integrate(1/self.__inv_eta**2*p0*p1)  
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
   
       def solve_A(self,u,p):  
          """  
          solves Av=f-Au-B^*p (v=0 on fixed_u_mask)  
          """  
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
          self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))  
          return  self.__pde_u.getSolution(verbose=self.show_details)  
   
   
       def solve_prec(self,p):  
         a=self.__inv_eta*p  
         return a-util.integrate(a)/self.vol  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
38    
39    class DarcyFlow(object):
40        """
41        Represents and solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44    
45        M{u_{i,i} = f}
46    
47        where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents
48        the permeability.
49    
50        @note: The problem is solved in a least squares formulation.
51        """
52    
53        def __init__(self, domain):
54            """
55            Initializes the Darcy flux problem.
56    
57            @param domain: domain of the problem
58            @type domain: L{Domain}
59            """
60            self.domain=domain
61            self.__pde_v=LinearPDESystem(domain)
62            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
63            self.__pde_v.setSymmetryOn()
64            self.__pde_p=LinearSinglePDE(domain)
65            self.__pde_p.setSymmetryOn()
66            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
67            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
68    
69        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70            """
71            Assigns values to model parameters.
72    
73            @param f: volumetric sources/sinks
74            @type f: scalar value on the domain, e.g. L{Data}
75            @param g: flux sources/sinks
76            @type g: vector value on the domain, e.g. L{Data}
77            @param location_of_fixed_pressure: mask for locations where pressure is fixed
78            @type location_of_fixed_pressure: scalar value on the domain, e.g. L{Data}
79            @param location_of_fixed_flux: mask for locations where flux is fixed
80            @type location_of_fixed_flux: vector value on the domain (e.g. L{Data})
81            @param permeability: permeability tensor. If scalar C{s} is given the
82                                 tensor with C{s} on the main diagonal is used. If
83                                 vector C{v} is given the tensor with C{v} on the
84                                 main diagonal is used.
85            @type permeability: scalar, vector or tensor values on the domain, e.g.
86                                L{Data}
87    
88            @note: the values of parameters which are not set by calling
89                   C{setValue} are not altered
90            @note: at any point on the boundary of the domain the pressure
91                   (C{location_of_fixed_pressure}) >0 or the normal component of
92                   the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
93                   the normal is along the M{x_i} axis.
94            """
95            if f !=None:
96               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
97               if f.isEmpty():
98                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
99               else:
100                   if f.getRank()>0: raise ValueError,"illegal rank of f."
101               self.f=f
102            if g !=None:
103               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
104               if g.isEmpty():
105                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
106               else:
107                 if not g.getShape()==(self.domain.getDim(),):
108                   raise ValueError,"illegal shape of g"
109               self.__g=g
110    
111            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
112            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
113    
114            if permeability!=None:
115               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
116               if perm.getRank()==0:
117                   perm=perm*util.kronecker(self.domain.getDim())
118               elif perm.getRank()==1:
119                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
120                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
121               elif perm.getRank()==2:
122                  pass
123               else:
124                  raise ValueError,"illegal rank of permeability."
125               self.__permeability=perm
126               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
127    
128    
129        def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
130            """
131            Returns the flux for a given pressure C{p}.
132    
133            The flux is equal to C{fixed_flux} on locations where
134            C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}
135            and C{f} are used.
136    
137            @param p: pressure
138            @type p: scalar value on the domain, e.g. L{Data}
139            @param fixed_flux: flux on the locations of the domain marked by
140                               C{location_of_fixed_flux}
141            @type fixed_flux: vector values on the domain, e.g. L{Data}
142            @param tol: relative tolerance to be used
143            @type tol: positive float
144            @return: flux
145            @rtype: L{Data}
146            @note: the method uses the least squares solution
147                   M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator
148                   and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
149            """
150            self.__pde_v.setTolerance(tol)
151            self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
152            return self.__pde_v.getSolution(verbose=show_details)
153    
154        def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
155            """
156            Solves the problem.
157    
158            The iteration is terminated if the error in the pressure is less than
159            M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
160            side (see escript user's guide for details).
161    
162            @param u0: initial guess for the flux. At locations in the domain
163                       marked by C{location_of_fixed_flux} the value of C{u0} is
164                       kept unchanged.
165            @type u0: vector value on the domain, e.g. L{Data}
166            @param p0: initial guess for the pressure. At locations in the domain
167                       marked by C{location_of_fixed_pressure} the value of C{p0}
168                       is kept unchanged.
169            @type p0: scalar value on the domain, e.g. L{Data}
170            @param atol: absolute tolerance for the pressure
171            @type atol: non-negative C{float}
172            @param rtol: relative tolerance for the pressure
173            @type rtol: non-negative C{float}
174            @param sub_rtol: tolerance to be used in the sub iteration. It is
175                             recommended that M{sub_rtol<rtol*5.e-3}
176            @type sub_rtol: positive-negative C{float}
177            @param verbose: if True information on iteration progress is printed
178            @type verbose: C{bool}
179            @param show_details: if True information on the sub-iteration process
180                                 is printed
181            @type show_details: C{bool}
182            @return: flux and pressure
183            @rtype: C{tuple} of L{Data}
184    
185            @note: the problem is solved in a least squares formulation:
186    
187            M{(I+D^*D)u+Qp=D^*f+g}
188    
189            M{Q^*u+Q^*Qp=Q^*g}
190    
191            where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the
192            permeability M{k_{ij}}. We eliminate the flux from the problem by
193            setting
194    
195            M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}
196    
197            from the first equation. Inserted into the second equation we get
198    
199            M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}
200            on C{location_of_fixed_pressure}
201    
202            which is solved using the PCG method (precondition is M{Q^*Q}).
203            In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}
204            need to be solved using a sub-iteration scheme.
205            """
206            self.verbose=verbose
207            self.show_details= show_details and self.verbose
208            self.__pde_v.setTolerance(sub_rtol)
209            self.__pde_p.setTolerance(sub_rtol)
210            u2=u0*self.__pde_v.getCoefficient("q")
211            #
212            # first the reference velocity is calculated from
213            #
214            #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
215            #
216            self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
217            u_ref=self.__pde_v.getSolution(verbose=show_details)
218            if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
219            self.__pde_v.setValue(r=Data())
220            #
221            #   and then we calculate a reference pressure
222            #
223            #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
224            #
225            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
226            p_ref=self.__pde_p.getSolution(verbose=self.show_details)
227            if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
228            self.__pde_p.setValue(r=Data())
229            #
230            #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
231            #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
232            #
233            #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
234            #
235            #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
236            #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
237            #
238            #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
239            #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
240            #
241            #  we use the (du,Qdp) to represent the resudual
242            #  Q^*Q is a preconditioner
243            #
244            #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
245            #
246            Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
247            norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
248            ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
249            if not ATOL>0:
250                raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)
251            if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
252            #
253            #   caclulate the initial residual
254            #
255            self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
256            du=self.__pde_v.getSolution(verbose=show_details)
257            r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
258            dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
259            util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
260            return u_ref+r[1],dp
261    
262        def __Aprod_PCG(self,p):
263            if self.show_details: print "DarcyFlux: Applying operator"
264            Qp=util.tensor_mult(self.__permeability,util.grad(p))
265            self.__pde_v.setValue(Y=Qp,X=Data())
266            w=self.__pde_v.getSolution(verbose=self.show_details)
267            return ArithmeticTuple(-Qp,w)
268    
269        def __inner_PCG(self,p,r):
270            a=util.tensor_mult(self.__permeability,util.grad(p))
271            out=-util.integrate(util.inner(a,r[0]+r[1]))
272            return out
273    
274        def __Msolve_PCG(self,r):
275            if self.show_details: print "DarcyFlux: Applying preconditioner"
276            self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))
277            return self.__pde_p.getSolution(verbose=self.show_details)
278    
279  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
280        """        """
281        solves        Represents and solves the problem
282    
283          M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}
284    
285            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0
                 u_{i,i}=0  
286    
287            u=0 where  fixed_u_mask>0        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}
           eta*(u_{i,j}+u_{j,i})*n_j=surface_stress  
288    
289        if surface_stress is not give 0 is assumed.        If C{surface_stress} is not given 0 is assumed.
290    
291        typical usage:        Typical usage::
292    
293              sp=StokesProblemCartesian(domain)            sp = StokesProblemCartesian(domain)
294              sp.setTolerance()            sp.setTolerance()
295              sp.initialize(...)            sp.initialize(...)
296              v,p=sp.solve(v0,p0)            v,p = sp.solve(v0,p0)
297        """        """
298        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
299             """
300             Initializes the Stokes Problem.
301    
302             @param domain: domain of the problem. The approximation order needs
303                            to be two.
304             @type domain: L{Domain}
305             @warning: The approximation order needs to be two otherwise you may
306                       see oscillations in the pressure.
307             """
308           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309           self.domain=domain           self.domain=domain
310           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
311           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())
312           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
313           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
314                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
315    
316           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
317           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
318             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
319           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
320    
321           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 174  class StokesProblemCartesian(Homogeneous Line 323  class StokesProblemCartesian(Homogeneous
323           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
324           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
325    
326        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(),stress=Data()):
327          self.eta=eta           """
328          A =self.__pde_u.createCoefficientOfGeneralPDE("A")           Assigns values to the model parameters.
329      self.__pde_u.setValue(A=Data())  
330          for i in range(self.domain.getDim()):           @param f: external force
331          for j in range(self.domain.getDim()):           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
332              A[i,j,j,i] += 1.           @param fixed_u_mask: mask of locations with fixed velocity
333              A[i,j,i,j] += 1.           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}
334      self.__pde_prec.setValue(D=1/self.eta)                               or similar
335          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)           @param eta: viscosity
336             @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar
337        def B(self,arg):           @param surface_stress: normal surface stress
338           d=util.div(arg)           @type surface_stress: L{Vector} object on L{FunctionSpace},
339           self.__pde_proj.setValue(Y=d)                                 L{FunctionOnBoundary} or similar
340             @param stress: initial stress
341             @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or
342                           similar
343             @note: All values need to be set.
344             """
345             self.eta=eta
346             A =self.__pde_u.createCoefficient("A")
347             self.__pde_u.setValue(A=Data())
348             for i in range(self.domain.getDim()):
349                 for j in range(self.domain.getDim()):
350                     A[i,j,j,i] += 1.
351                     A[i,j,i,j] += 1.
352             self.__pde_prec.setValue(D=1/self.eta)
353             self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
354             self.__stress=stress
355    
356          def B(self,v):
357             """
358             Returns M{div(v)}.
359             @return: M{div(v)}
360             @rtype: equal to the type of p
361    
362             @note: Boundary conditions on p should be zero!
363             """
364             if self.show_details: print "apply divergence:"
365             self.__pde_proj.setValue(Y=-util.div(v))
366           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setTolerance(self.getSubProblemTolerance())
367           return self.__pde_proj.getSolution(verbose=self.show_details)           return self.__pde_proj.getSolution(verbose=self.show_details)
368    
369        def inner(self,p0,p1):        def inner_pBv(self,p,Bv):
370           s0=util.interpolate(p0,Function(self.domain))           """
371           s1=util.interpolate(p1,Function(self.domain))           Returns inner product of element p and Bv (overwrite).
372    
373             @type p: equal to the type of p
374             @type Bv: equal to the type of result of operator B
375             @return: inner product of p and Bv
376             @rtype: equal to the type of p
377             """
378             s0=util.interpolate(p,Function(self.domain))
379             s1=util.interpolate(Bv,Function(self.domain))
380           return util.integrate(s0*s1)           return util.integrate(s0*s1)
381    
382        def inner_a(self,a0,a1):        def inner_p(self,p0,p1):
383           p0=util.interpolate(a0[1],Function(self.domain))           """
384           p1=util.interpolate(a1[1],Function(self.domain))           Returns inner product of element p0 and p1 (overwrite).
          alfa=(1/self.vol)*util.integrate(p0)  
          beta=(1/self.vol)*util.integrate(p1)  
      v0=util.grad(a0[0])  
      v1=util.grad(a1[0])  
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
385    
386             @type p0: equal to the type of p
387             @type p1: equal to the type of p
388             @return: inner product of p0 and p1
389             @rtype: equal to the type of p
390             """
391             s0=util.interpolate(p0/self.eta,Function(self.domain))
392             s1=util.interpolate(p1/self.eta,Function(self.domain))
393             return util.integrate(s0*s1)
394    
395          def inner_v(self,v0,v1):
396             """
397             Returns inner product of two elements v0 and v1 (overwrite).
398    
399        def getStress(self,u):           @type v0: equal to the type of v
400           mg=util.grad(u)           @type v1: equal to the type of v
401           return 2.*self.eta*util.symmetric(mg)           @return: inner product of v0 and v1
402        def getEtaEffective(self):           @rtype: equal to the type of v
403           return self.eta           """
404             gv0=util.grad(v0)
405             gv1=util.grad(v1)
406             return util.integrate(util.inner(gv0,gv1))
407    
408        def solve_A(self,u,p):        def solve_A(self,u,p):
409           """           """
410           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).
411           """           """
412             if self.show_details: print "solve for velocity:"
413           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
414           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           if self.__stress.isEmpty():
415           return  self.__pde_u.getSolution(verbose=self.show_details)              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
416             else:
417                self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
418             out=self.__pde_u.getSolution(verbose=self.show_details)
419             return out
420    
421        def solve_prec(self,p):        def solve_prec(self,p):
422       #proj=Projector(domain=self.domain, reduce = True, fast=False)           """
423           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           Applies the preconditioner.
424           self.__pde_prec.setValue(Y=p)           """
425           q=self.__pde_prec.getSolution(verbose=self.show_details)           if self.show_details: print "apply preconditioner:"
          return q  
   
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
426           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
427           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
428           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
      q0=util.interpolate(q,Function(self.domain))  
          print util.inf(q*q0),util.sup(q*q0)  
          q-=(1/self.vol)*util.integrate(q0)  
          print util.inf(q*q0),util.sup(q*q0)  
429           return q           return q
430    
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
   
   

Legend:
Removed from v.1673  
changed lines
  Added in v.2169

  ViewVC Help
Powered by ViewVC 1.1.26