/[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 2474 by gross, Tue Jun 16 06:32:15 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, SolverOptions
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, adaptSubTolerance=True):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56        @param useReduced: uses reduced oreder on flux and pressure
57        @type useReduced: C{bool}
58        @param adaptSubTolerance: switches on automatic subtolerance selection
59        @type adaptSubTolerance: C{bool}    
60            """
61            self.domain=domain
62            if weight == None:
63               s=self.domain.getSize()
64               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65            else:
66               self.__l=weight
67            self.__pde_v=LinearPDESystem(domain)
68            if useReduced: self.__pde_v.setReducedOrderOn()
69            self.__pde_v.setSymmetryOn()
70            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
71            self.__pde_p=LinearSinglePDE(domain)
72            self.__pde_p.setSymmetryOn()
73            if useReduced: self.__pde_p.setReducedOrderOn()
74            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
75            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
76            self.setTolerance()
77            self.setAbsoluteTolerance()
78        self.__adaptSubTolerance=adaptSubTolerance
79        self.verbose=False
80        def getSolverOptionsFlux(self):
81        """
82        Returns the solver options used to solve the flux problems
83        
84        M{(I+D^*D)u=F}
85        
86        @return: L{SolverOptions}
87        """
88        return self.__pde_v.getSolverOptions()
89        def setSolverOptionsFlux(self, options=None):
90        """
91        Sets the solver options used to solve the flux problems
92        
93        M{(I+D^*D)u=F}
94        
95        If C{options} is not present, the options are reset to default
96        @param options: L{SolverOptions}
97        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
98        """
99        return self.__pde_v.setSolverOptions(options)
100        def getSolverOptionsPressure(self):
101        """
102        Returns the solver options used to solve the pressure problems
103        
104        M{(Q^*Q)p=Q^*G}
105        
106        @return: L{SolverOptions}
107        """
108        return self.__pde_p.getSolverOptions()
109        def setSolverOptionsPressure(self, options=None):
110        """
111        Sets the solver options used to solve the pressure problems
112        
113        M{(Q^*Q)p=Q^*G}
114        
115        If C{options} is not present, the options are reset to default
116        @param options: L{SolverOptions}
117        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
118        """
119        return self.__pde_p.setSolverOptions(options)
120    
121        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122            """
123            assigns values to model parameters
124    
125            @param f: volumetic sources/sinks
126            @type f: scalar value on the domain (e.g. L{Data})
127            @param g: flux sources/sinks
128            @type g: vector values on the domain (e.g. L{Data})
129            @param location_of_fixed_pressure: mask for locations where pressure is fixed
130            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
131            @param location_of_fixed_flux:  mask for locations where flux is fixed.
132            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
133            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
134                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
135                                 C{v} on the main diagonal is used.
136            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
137    
138            @note: the values of parameters which are not set by calling C{setValue} are not altered.
139            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
140                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
141                   is along the M{x_i} axis.
142            """
143            if f !=None:
144               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
145               if f.isEmpty():
146                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
147               else:
148                   if f.getRank()>0: raise ValueError,"illegal rank of f."
149               self.__f=f
150            if g !=None:
151               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
152               if g.isEmpty():
153                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
154               else:
155                 if not g.getShape()==(self.domain.getDim(),):
156                   raise ValueError,"illegal shape of g"
157               self.__g=g
158    
159            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
160            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
161    
162            if permeability!=None:
163               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
164               if perm.getRank()==0:
165                   perm=perm*util.kronecker(self.domain.getDim())
166               elif perm.getRank()==1:
167                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
168                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
169               elif perm.getRank()==2:
170                  pass
171               else:
172                  raise ValueError,"illegal rank of permeability."
173               self.__permeability=perm
174               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
175    
176        def setTolerance(self,rtol=1e-4):
177            """
178            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
179    
180            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
181    
182            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}}.
183    
184            @param rtol: relative tolerance for the pressure
185            @type rtol: non-negative C{float}
186            """
187            if rtol<0:
188                raise ValueError,"Relative tolerance needs to be non-negative."
189            self.__rtol=rtol
190        def getTolerance(self):
191            """
192            returns the relative tolerance
193    
194            @return: current relative tolerance
195            @rtype: C{float}
196            """
197            return self.__rtol
198    
199        def setAbsoluteTolerance(self,atol=0.):
200            """
201            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
202    
203            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
204    
205            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}}.
206    
207            @param atol: absolute tolerance for the pressure
208            @type atol: non-negative C{float}
209            """
210            if atol<0:
211                raise ValueError,"Absolute tolerance needs to be non-negative."
212            self.__atol=atol
213        def getAbsoluteTolerance(self):
214           """
215           returns the absolute tolerance
216          
217           @return: current absolute tolerance
218           @rtype: C{float}
219           """
220           return self.__atol
221        def getSubProblemTolerance(self):
222        """
223        Returns a suitable subtolerance
224        @type: C{float}
225        """
226        return max(util.EPSILON**(0.75),self.getTolerance()**2)
227        def setSubProblemTolerance(self):
228             """
229             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
230             """
231         if self.__adaptSubTolerance:
232             sub_tol=self.getSubProblemTolerance()
233                 self.getSolverOptionsFlux().setTolerance(sub_tol)
234             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
235             self.getSolverOptionsPressure().setTolerance(sub_tol)
236             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
237             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
238    
239            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
240                  u_{i,i}=0           """
241             solves the problem.
242    
243            u=0 where  fixed_u_mask>0           The iteration is terminated if the residual norm is less then self.getTolerance().
           eta*(u_{i,j}+u_{j,i})*n_j=surface_stress  
244    
245        if surface_stress is not give 0 is assumed.           @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.
246             @type u0: vector value on the domain (e.g. L{Data}).
247             @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.
248             @type p0: scalar value on the domain (e.g. L{Data}).
249             @param verbose: if set some information on iteration progress are printed
250             @type verbose: C{bool}
251             @return: flux and pressure
252             @rtype: C{tuple} of L{Data}.
253    
254        typical usage:           @note: The problem is solved as a least squares form
255    
256              sp=StokesProblemCartesian(domain)           M{(I+D^*D)u+Qp=D^*f+g}
257              sp.setTolerance()           M{Q^*u+Q^*Qp=Q^*g}
             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)  
258    
259        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
260          self.eta=eta           We eliminate the flux form the problem by setting
261          A =self.__pde_u.createCoefficientOfGeneralPDE("A")  
262      self.__pde_u.setValue(A=Data())           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
         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)  
263    
264          # self.__pde_proj.setValue(D=1/eta)           form the first equation. Inserted into the second equation we get
         # 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  
265    
266             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
267    
268             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
269             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
270             """
271             self.verbose=verbose
272             rtol=self.getTolerance()
273             atol=self.getAbsoluteTolerance()
274         self.setSubProblemTolerance()
275        
276             num_corrections=0
277             converged=False
278             p=p0
279             norm_r=None
280             while not converged:
281                   v=self.getFlux(p, fixed_flux=u0)
282                   Qp=self.__Q(p)
283                   norm_v=self.__L2(v)
284                   norm_Qp=self.__L2(Qp)
285                   if norm_v == 0.:
286                      if norm_Qp == 0.:
287                         return v,p
288                      else:
289                        fac=norm_Qp
290                   else:
291                      if norm_Qp == 0.:
292                        fac=norm_v
293                      else:
294                        fac=2./(1./norm_v+1./norm_Qp)
295                   ATOL=(atol+rtol*fac)
296                   if self.verbose:
297                        print "DarcyFlux: L2 norm of v = %e."%norm_v
298                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
299                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
300                   if norm_r == None or norm_r>ATOL:
301                       if num_corrections>max_num_corrections:
302                             raise ValueError,"maximum number of correction steps reached."
303                       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)
304                       num_corrections+=1
305                   else:
306                       converged=True
307             return v,p
308        def __L2(self,v):
309             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
310    
311        def __Q(self,p):
312              return util.tensor_mult(self.__permeability,util.grad(p))
313    
314        def __Aprod(self,dp):
315              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
316              Qdp=self.__Q(dp)
317              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
318              du=self.__pde_v.getSolution()
319              # self.__pde_v.getOperator().saveMM("proj.mm")
320              return Qdp+du
321        def __inner_GMRES(self,r,s):
322             return util.integrate(util.inner(r,s))
323    
324        def __inner_PCG(self,p,r):
325             return util.integrate(util.inner(self.__Q(p), r))
326    
327        def __Msolve_PCG(self,r):
328          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
329              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
330              # self.__pde_p.getOperator().saveMM("prec.mm")
331              return self.__pde_p.getSolution()
332    
333        def getFlux(self,p=None, fixed_flux=Data()):
334            """
335            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
336            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
337            Note that C{g} and C{f} are used, see L{setValue}.
338    
339            @param p: pressure.
340            @type p: scalar value on the domain (e.g. L{Data}).
341            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
342            @type fixed_flux: vector values on the domain (e.g. L{Data}).
343            @param tol: relative tolerance to be used.
344            @type tol: positive C{float}.
345            @return: flux
346            @rtype: L{Data}
347            @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}}
348                   for the permeability M{k_{ij}}
349            """
350        self.setSubProblemTolerance()
351            g=self.__g
352            f=self.__f
353            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
354            if p == None:
355               self.__pde_v.setValue(Y=g)
356            else:
357               self.__pde_v.setValue(Y=g-self.__Q(p))
358            return self.__pde_v.getSolution()
359    
360  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
361        """       """
362        solves       solves
363    
364            -(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}
365                  u_{i,i}=0                  u_{i,i}=0
366    
367            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
368            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
369    
370        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
371    
372        typical usage:       typical usage:
373    
374              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
375              sp.setTolerance()              sp.setTolerance()
376              sp.initialize(...)              sp.initialize(...)
377              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
378        """       """
379        def __init__(self,domain,**kwargs):       def __init__(self,domain,adaptSubTolerance=True, **kwargs):
380           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           """
381             initialize the Stokes Problem
382    
383             @param domain: domain of the problem. The approximation order needs to be two.
384             @type domain: L{Domain}
385         @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
386         @type adaptSubTolerance: C{bool}
387             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
388             """
389             HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
390           self.domain=domain           self.domain=domain
391           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
392           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())
393           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
394           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)      
               
395           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
396           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
397           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
398    
399           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
400           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
401         self.__pde_proj.setValue(D=1)
402           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
403    
404        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
405             """
406         returns the solver options used  solve the equation for velocity.
407        
408         @rtype: L{SolverOptions}
409         """
410         return self.__pde_u.getSolverOptions()
411         def setSolverOptionsVelocity(self, options=None):
412             """
413         set the solver options for solving the equation for velocity.
414        
415         @param options: new solver  options
416         @type options: L{SolverOptions}
417         """
418             self.__pde_u.setSolverOptions(options)
419         def getSolverOptionsPressure(self):
420             """
421         returns the solver options used  solve the equation for pressure.
422         @rtype: L{SolverOptions}
423         """
424         return self.__pde_prec.getSolverOptions()
425         def setSolverOptionsPressure(self, options=None):
426             """
427         set the solver options for solving the equation for pressure.
428         @param options: new solver  options
429         @type options: L{SolverOptions}
430         """
431         self.__pde_prec.setSolverOptions(options)
432    
433         def setSolverOptionsDiv(self, options=None):
434             """
435         set the solver options for solving the equation to project the divergence of
436         the velocity onto the function space of presure.
437        
438         @param options: new solver options
439         @type options: L{SolverOptions}
440         """
441         self.__pde_prec.setSolverOptions(options)
442         def getSolverOptionsDiv(self):
443             """
444         returns the solver options for solving the equation to project the divergence of
445         the velocity onto the function space of presure.
446        
447         @rtype: L{SolverOptions}
448         """
449         return self.__pde_prec.getSolverOptions()
450         def setSubProblemTolerance(self):
451             """
452         Updates the tolerance for subproblems
453             """
454         if self.adaptSubTolerance():
455                 sub_tol=self.getSubProblemTolerance()
456             self.getSolverOptionsDiv().setTolerance(sub_tol)
457             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
458             self.getSolverOptionsPressure().setTolerance(sub_tol)
459             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
460             self.getSolverOptionsVelocity().setTolerance(sub_tol)
461             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
462            
463    
464         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
465            """
466            assigns values to the model parameters
467    
468            @param f: external force
469            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
470            @param fixed_u_mask: mask of locations with fixed velocity.
471            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
472            @param eta: viscosity
473            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
474            @param surface_stress: normal surface stress
475            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
476            @param stress: initial stress
477        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
478            @note: All values needs to be set.
479            """
480          self.eta=eta          self.eta=eta
481          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
482      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
483          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
484          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
485              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
486              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
487      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
488          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)
489            self.__f=f
490            self.__surface_stress=surface_stress
491            self.__stress=stress
492    
493         def Bv(self,v):
494             """
495             returns inner product of element p and div(v)
496    
497             @param p: a pressure increment
498             @param v: a residual
499             @return: inner product of element p and div(v)
500             @rtype: C{float}
501             """
502             self.__pde_proj.setValue(Y=-util.div(v))
503             return self.__pde_proj.getSolution()
504    
505         def inner_pBv(self,p,Bv):
506             """
507             returns inner product of element p and Bv=-div(v)
508    
509             @param p: a pressure increment
510             @param v: a residual
511             @return: inner product of element p and Bv=-div(v)
512             @rtype: C{float}
513             """
514             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
515    
516        def B(self,arg):       def inner_p(self,p0,p1):
517           d=util.div(arg)           """
518           self.__pde_proj.setValue(Y=d)           Returns inner product of p0 and p1
519           self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
520           return self.__pde_proj.getSolution(verbose=self.show_details)           @param p0: a pressure
521             @param p1: a pressure
522        def inner(self,p0,p1):           @return: inner product of p0 and p1
523           s0=util.interpolate(p0,Function(self.domain))           @rtype: C{float}
524           s1=util.interpolate(p1,Function(self.domain))           """
525             s0=util.interpolate(p0/self.eta,Function(self.domain))
526             s1=util.interpolate(p1/self.eta,Function(self.domain))
527           return util.integrate(s0*s1)           return util.integrate(s0*s1)
528    
529        def inner_a(self,a0,a1):       def norm_v(self,v):
530           p0=util.interpolate(a0[1],Function(self.domain))           """
531           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
532           alfa=(1/self.vol)*util.integrate(p0)  
533           beta=(1/self.vol)*util.integrate(p1)           @param v: a velovity
534       v0=util.grad(a0[0])           @return: norm of v
535       v1=util.grad(a1[0])           @rtype: non-negative C{float}
536           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           """
537             return util.sqrt(util.integrate(util.length(util.grad(v))))
538    
539        def getStress(self,u):       def getV(self, p, v0):
540           mg=util.grad(u)           """
541           return 2.*self.eta*util.symmetric(mg)           return the value for v for a given p (overwrite)
       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)-p*util.kronecker(self.domain))  
          return  self.__pde_u.getSolution(verbose=self.show_details)  
   
   
       def solve_prec(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)  
          return q  
   
       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  
542    
543             @param p: a pressure
544             @param v0: a initial guess for the value v to return.
545             @return: v given as M{v= A^{-1} (f-B^*p)}
546             """
547             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
548             if self.__stress.isEmpty():
549                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
550             else:
551                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
552             out=self.__pde_u.getSolution()
553             return  out
554    
555         def norm_Bv(self,Bv):
556            """
557            Returns Bv (overwrite).
558    
559            @rtype: equal to the type of p
560            @note: boundary conditions on p should be zero!
561            """
562            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
563    
564         def solve_AinvBt(self,p):
565             """
566             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
567    
568             @param p: a pressure increment
569             @return: the solution of M{Av=B^*p}
570             @note: boundary conditions on v should be zero!
571             """
572             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
573             out=self.__pde_u.getSolution()
574             return  out
575    
576         def solve_prec(self,Bv):
577             """
578             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
579             with accuracy L{self.getSubProblemTolerance()}
580    
581             @param v: velocity increment
582             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
583             @note: boundary conditions on p are zero.
584             """
585             self.__pde_prec.setValue(Y=Bv)
586             return self.__pde_prec.getSolution()

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

  ViewVC Help
Powered by ViewVC 1.1.26