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

Diff of /trunk/escript/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 2386 by gross, Wed Apr 15 03:54:25 2009 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__="https://launchpad.net/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, GMRES
38    
39  class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  class DarcyFlow(object):
40        """      """
41        solves      solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        M{u_{i,i} = f}
45    
46        where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47    
48        @note: The problem is solved in a least squares formulation.
49        """
50    
51        def __init__(self, domain, weight=None, useReduced=False):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56            """
57            self.domain=domain
58            if weight == None:
59               s=self.domain.getSize()
60               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61            else:
62               self.__l=weight
63            self.__pde_v=LinearPDESystem(domain)
64            if useReduced: self.__pde_v.setReducedOrderOn()
65            self.__pde_v.setSymmetryOn()
66            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68            self.__pde_p=LinearSinglePDE(domain)
69            self.__pde_p.setSymmetryOn()
70            if useReduced: self.__pde_p.setReducedOrderOn()
71            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73            self.setTolerance()
74            self.setAbsoluteTolerance()
75            self.setSubProblemTolerance()
76    
77        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
78            """
79            assigns values to model parameters
80    
81            @param f: volumetic sources/sinks
82            @type f: scalar value on the domain (e.g. L{Data})
83            @param g: flux sources/sinks
84            @type g: vector values on the domain (e.g. L{Data})
85            @param location_of_fixed_pressure: mask for locations where pressure is fixed
86            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
87            @param location_of_fixed_flux:  mask for locations where flux is fixed.
88            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
89            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91                                 C{v} on the main diagonal is used.
92            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93    
94            @note: the values of parameters which are not set by calling C{setValue} are not altered.
95            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
97                   is along the M{x_i} axis.
98            """
99            if f !=None:
100               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101               if f.isEmpty():
102                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103               else:
104                   if f.getRank()>0: raise ValueError,"illegal rank of f."
105               self.__f=f
106            if g !=None:
107               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108               if g.isEmpty():
109                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110               else:
111                 if not g.getShape()==(self.domain.getDim(),):
112                   raise ValueError,"illegal shape of g"
113               self.__g=g
114    
115            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117    
118            if permeability!=None:
119               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120               if perm.getRank()==0:
121                   perm=perm*util.kronecker(self.domain.getDim())
122               elif perm.getRank()==1:
123                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125               elif perm.getRank()==2:
126                  pass
127               else:
128                  raise ValueError,"illegal rank of permeability."
129               self.__permeability=perm
130               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
131    
132        def setTolerance(self,rtol=1e-4):
133            """
134            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135    
136            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137    
138            where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
139    
140            @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142            """
143            if rtol<0:
144                raise ValueError,"Relative tolerance needs to be non-negative."
145            self.__rtol=rtol
146        def getTolerance(self):
147            """
148            returns the relative tolerance
149    
150            @return: current relative tolerance
151            @rtype: C{float}
152            """
153            return self.__rtol
154    
155        def setAbsoluteTolerance(self,atol=0.):
156            """
157            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158    
159            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160    
161            where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
162    
163            @param atol: absolute tolerance for the pressure
164            @type atol: non-negative C{float}
165            """
166            if atol<0:
167                raise ValueError,"Absolute tolerance needs to be non-negative."
168            self.__atol=atol
169        def getAbsoluteTolerance(self):
170           """
171           returns the absolute tolerance
172          
173           @return: current absolute tolerance
174           @rtype: C{float}
175           """
176           return self.__atol
177    
178            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i      def setSubProblemTolerance(self,rtol=None):
179                  u_{i,i}=0           """
180             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181             C{self.getTolerance()**2} is used.
182    
183            u=0 where  fixed_u_mask>0           @param rtol: relative tolerence
184            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress           @type rtol: positive C{float}
185             """
186             if rtol == None:
187                  if self.getTolerance()<=0.:
188                      raise ValueError,"A positive relative tolerance must be set."
189                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190             else:
191                 if rtol<=0:
192                     raise ValueError,"sub-problem tolerance must be positive."
193                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194    
195        if surface_stress is not give 0 is assumed.      def getSubProblemTolerance(self):
196             """
197             Returns the subproblem reduction factor.
198    
199        typical usage:           @return: subproblem reduction factor
200             @rtype: C{float}
201             """
202             return self.__sub_tol
203    
204              sp=StokesProblemCartesian(domain)      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205              sp.setTolerance()           """
206              sp.initialize(...)           solves the problem.
             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)  
207    
208        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):           The iteration is terminated if the residual norm is less then self.getTolerance().
         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)  
209    
210          # self.__pde_proj.setValue(D=1/eta)           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
211          # self.__pde_proj.setValue(Y=1.)           @type u0: vector value on the domain (e.g. L{Data}).
212          # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
213          self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))           @type p0: scalar value on the domain (e.g. L{Data}).
214             @param verbose: if set some information on iteration progress are printed
215        def B(self,arg):           @type verbose: C{bool}
216           a=util.div(arg, ReducedFunction(self.domain))           @param show_details:  if set information on the subiteration process are printed.
217           return a-util.integrate(a)/self.vol           @type show_details: C{bool}
218             @return: flux and pressure
219        def inner(self,p0,p1):           @rtype: C{tuple} of L{Data}.
          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  
220    
221        def solve_A(self,u,p):           @note: The problem is solved as a least squares form
222           """  
223           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           M{(I+D^*D)u+Qp=D^*f+g}
224           """           M{Q^*u+Q^*Qp=Q^*g}
          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)  
225    
226             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227             We eliminate the flux form the problem by setting
228    
229        def solve_prec(self,p):           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
         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  
230    
231             form the first equation. Inserted into the second equation we get
232    
233             M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure
234    
235             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
236             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237             """
238             self.verbose=verbose
239             self.show_details= show_details and self.verbose
240             rtol=self.getTolerance()
241             atol=self.getAbsoluteTolerance()
242             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243    
244             num_corrections=0
245             converged=False
246             p=p0
247             norm_r=None
248             while not converged:
249                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250                   Qp=self.__Q(p)
251                   norm_v=self.__L2(v)
252                   norm_Qp=self.__L2(Qp)
253                   if norm_v == 0.:
254                      if norm_Qp == 0.:
255                         return v,p
256                      else:
257                        fac=norm_Qp
258                   else:
259                      if norm_Qp == 0.:
260                        fac=norm_v
261                      else:
262                        fac=2./(1./norm_v+1./norm_Qp)
263                   ATOL=(atol+rtol*fac)
264                   if self.verbose:
265                        print "DarcyFlux: L2 norm of v = %e."%norm_v
266                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268                   if norm_r == None or norm_r>ATOL:
269                       if num_corrections>max_num_corrections:
270                             raise ValueError,"maximum number of correction steps reached."
271                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272                       num_corrections+=1
273                   else:
274                       converged=True
275             return v,p
276        def __L2(self,v):
277             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278    
279        def __Q(self,p):
280              return util.tensor_mult(self.__permeability,util.grad(p))
281    
282        def __Aprod(self,dp):
283              self.__pde_v.setTolerance(self.getSubProblemTolerance())
284              if self.show_details: print "DarcyFlux: Applying operator"
285              Qdp=self.__Q(dp)
286              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288              # self.__pde_v.getOperator().saveMM("proj.mm")
289              return Qdp+du
290        def __inner_GMRES(self,r,s):
291             return util.integrate(util.inner(r,s))
292    
293        def __inner_PCG(self,p,r):
294             return util.integrate(util.inner(self.__Q(p), r))
295    
296        def __Msolve_PCG(self,r):
297              self.__pde_p.setTolerance(self.getSubProblemTolerance())
298              if self.show_details: print "DarcyFlux: Applying preconditioner"
299              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300              # self.__pde_p.getOperator().saveMM("prec.mm")
301              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302    
303        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304            """
305            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
306            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
307            Note that C{g} and C{f} are used, see L{setValue}.
308    
309            @param p: pressure.
310            @type p: scalar value on the domain (e.g. L{Data}).
311            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
312            @type fixed_flux: vector values on the domain (e.g. L{Data}).
313            @param tol: relative tolerance to be used.
314            @type tol: positive C{float}.
315            @return: flux
316            @rtype: L{Data}
317            @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
318                   for the permeability M{k_{ij}}
319            """
320            self.__pde_v.setTolerance(self.getSubProblemTolerance())
321            g=self.__g
322            f=self.__f
323            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
324            if p == None:
325               self.__pde_v.setValue(Y=g)
326            else:
327               self.__pde_v.setValue(Y=g-self.__Q(p))
328            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329    
330  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
331        """       """
332        solves       solves
333    
334            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
335                  u_{i,i}=0                  u_{i,i}=0
336    
337            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
338            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
339    
340        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
341    
342        typical usage:       typical usage:
343    
344              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
345              sp.setTolerance()              sp.setTolerance()
346              sp.initialize(...)              sp.initialize(...)
347              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
348        """       """
349        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
350             """
351             initialize the Stokes Problem
352    
353             @param domain: domain of the problem. The approximation order needs to be two.
354             @type domain: L{Domain}
355             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
356             """
357           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
358           self.domain=domain           self.domain=domain
359           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
360           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())
361           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
362             # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364                
365           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
366           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
367             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
369    
370           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
371           self.__pde_proj.setReducedOrderOn()          """
372           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
373           self.__pde_proj.setValue(D=1.)  
374            @param f: external force
375            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
376            @param fixed_u_mask: mask of locations with fixed velocity.
377            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
378            @param eta: viscosity
379            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
380            @param surface_stress: normal surface stress
381            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
382            @param stress: initial stress
383        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
384            @note: All values needs to be set.
385    
386        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          """
387          self.eta=eta          self.eta=eta
388          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
389      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
390          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
391          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
392              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
393              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
394      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
395          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)
396            self.__f=f
397            self.__surface_stress=surface_stress
398            self.__stress=stress
399    
400        def B(self,arg):       def inner_pBv(self,p,v):
401           d=util.div(arg)           """
402           self.__pde_proj.setValue(Y=d)           returns inner product of element p and div(v)
403           self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
404           return self.__pde_proj.getSolution(verbose=self.show_details)           @param p: a pressure increment
405             @param v: a residual
406        def inner(self,p0,p1):           @return: inner product of element p and div(v)
407           s0=util.interpolate(p0,Function(self.domain))           @rtype: C{float}
408           s1=util.interpolate(p1,Function(self.domain))           """
409             return util.integrate(-p*util.div(v))
410    
411         def inner_p(self,p0,p1):
412             """
413             Returns inner product of p0 and p1
414    
415             @param p0: a pressure
416             @param p1: a pressure
417             @return: inner product of p0 and p1
418             @rtype: C{float}
419             """
420             s0=util.interpolate(p0/self.eta,Function(self.domain))
421             s1=util.interpolate(p1/self.eta,Function(self.domain))
422           return util.integrate(s0*s1)           return util.integrate(s0*s1)
423    
424        def inner_a(self,a0,a1):       def norm_v(self,v):
425           p0=util.interpolate(a0[1],Function(self.domain))           """
426           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
427           alfa=(1/self.vol)*util.integrate(p0)  
428           beta=(1/self.vol)*util.integrate(p1)           @param v: a velovity
429       v0=util.grad(a0[0])           @return: norm of v
430       v1=util.grad(a1[0])           @rtype: non-negative C{float}
431           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           """
432             return util.sqrt(util.integrate(util.length(util.grad(v))))
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
433    
434        def solve_A(self,u,p):       def getV(self, p, v0):
435           """           """
436           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
437    
438             @param p: a pressure
439             @param v0: a initial guess for the value v to return.
440             @return: v given as M{v= A^{-1} (f-B^*p)}
441           """           """
442           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
443           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
444           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
445                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
446             else:
447                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
448             out=self.__pde_u.getSolution(verbose=self.show_details)
449             return  out
450    
451    
452        def solve_prec(self,p):           raise NotImplementedError,"no v calculation implemented."
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
453    
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          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)  
          return q  
   
       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  
454    
455         def norm_Bv(self,v):
456            """
457            Returns Bv (overwrite).
458    
459            @rtype: equal to the type of p
460            @note: boundary conditions on p should be zero!
461            """
462            return util.sqrt(util.integrate(util.div(v)**2))
463    
464         def solve_AinvBt(self,p):
465             """
466             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
467    
468             @param p: a pressure increment
469             @return: the solution of M{Av=B^*p}
470             @note: boundary conditions on v should be zero!
471             """
472             self.__pde_u.setTolerance(self.getSubProblemTolerance())
473             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
474             out=self.__pde_u.getSolution(verbose=self.show_details)
475             return  out
476    
477         def solve_precB(self,v):
478             """
479             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
480             with accuracy L{self.getSubProblemTolerance()} (overwrite).
481    
482             @param v: velocity increment
483             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
484             @note: boundary conditions on p are zero.
485             """
486             self.__pde_prec.setValue(Y=-util.div(v))
487             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
488             return self.__pde_prec.getSolution(verbose=self.show_details)

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

  ViewVC Help
Powered by ViewVC 1.1.26