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

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

  ViewVC Help
Powered by ViewVC 1.1.26