/[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 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 3071 by gross, Wed Jul 21 05:37:30 2010 UTC
# Line 1  Line 1 
1  # $Id:$  # -*- coding: utf-8 -*-
2    ########################################################
3  #  #
4  #######################################################  # Copyright (c) 2003-2010 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7  #  #
8  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
9  #  # Licensed under the Open Software License version 3.0
10  #                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  
 #  
 #######################################################  
11  #  #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
34    
35  from escript import *  from escript import *
36  import util  import util
37  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class DarcyFlow(object):
41       """
42       solves the problem
43      
44       *u_i+k_{ij}*p_{,j} = g_i*
45       *u_{i,i} = f*
46      
47       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49       :note: The problem is solved in a least squares formulation.
50       """
51      
52       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):
53          """
54          initializes the Darcy flux problem
55          :param domain: domain of the problem
56          :type domain: `Domain`
57          :param useReduced: uses reduced oreder on flux and pressure
58          :type useReduced: ``bool``
59          :param adaptSubTolerance: switches on automatic subtolerance selection
60          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          """
64          self.domain=domain
65          self.solveForFlux=solveForFlux
66          self.useReduced=useReduced
67          self.__adaptSubTolerance=adaptSubTolerance
68          self.verbose=False
69          
70          self.__pde_k=LinearPDESystem(domain)
71          self.__pde_k.setSymmetryOn()
72          if self.useReduced: self.__pde_k.setReducedOrderOn()
73    
74          self.__pde_p=LinearSinglePDE(domain)
75          self.__pde_p.setSymmetryOn()
76          if self.useReduced: self.__pde_p.setReducedOrderOn()
77    
78          self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting
79          # self.__pde_l.setSymmetryOn()
80          # if self.useReduced: self.__pde_l.setReducedOrderOn()
81          self.setTolerance()
82          self.setAbsoluteTolerance()
83          self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84          self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85          
86       def getSolverOptionsFlux(self):
87          """
88          Returns the solver options used to solve the flux problems
89          
90          *K^{-1} u=F*
91          
92          :return: `SolverOptions`
93          """
94          return self.__pde_k.getSolverOptions()
95          
96       def setSolverOptionsFlux(self, options=None):
97          """
98          Sets the solver options used to solve the flux problems
99          
100          *K^{-1}u=F*
101          
102          If ``options`` is not present, the options are reset to default
103          
104          :param options: `SolverOptions`
105          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
106          """
107          return self.__pde_v.setSolverOptions(options)
108        
109       def getSolverOptionsPressure(self):
110          """
111          Returns the solver options used to solve the pressure problems
112          
113          *(Q^* K Q)p=-Q^*G*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_p.getSolverOptions()
118          
119       def setSolverOptionsPressure(self, options=None):
120        """        """
121        solves        Sets the solver options used to solve the pressure problems
122          
123          *(Q^* K Q)p=-Q^*G*
124          
125          If ``options`` is not present, the options are reset to default
126          
127          :param options: `SolverOptions`
128          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
129          """
130          return self.__pde_p.setSolverOptions(options)
131    
132       def getSolverOptionsWeighting(self):
133          """
134          Returns the solver options used to solve the pressure problems
135    
136          *(D K D^*)p=-D F*
137    
138          :return: `SolverOptions`
139          """
140          return self.__pde_l.getSolverOptions()
141    
142       def setSolverOptionsWeighting(self, options=None):
143          """
144          Sets the solver options used to solve the pressure problems
145    
146          *(D K D^*)p=-D F*
147    
148          If ``options`` is not present, the options are reset to default
149    
150          :param options: `SolverOptions`
151          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
152          """
153          return self.__pde_l.setSolverOptions(options)
154    
155    
156       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
157          """
158          assigns values to model parameters
159          
160          :param f: volumetic sources/sinks
161          :type f: scalar value on the domain (e.g. `Data`)
162          :param g: flux sources/sinks
163          :type g: vector values on the domain (e.g. `Data`)
164          :param location_of_fixed_pressure: mask for locations where pressure is fixed
165          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
166          :param location_of_fixed_flux:  mask for locations where flux is fixed.
167          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
168          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
169          ``s`` on the main diagonal is used.
170          :type permeability: scalar or tensor values on the domain (e.g. `Data`)
171          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
172          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
173          or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
174          is along the *x_i* axis.
175          """
176          if f !=None:
177         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178         if f.isEmpty():
179            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180         else:
181            if f.getRank()>0: raise ValueError,"illegal rank of f."
182            self.__f=f
183          if g !=None:
184         g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185         if g.isEmpty():
186            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         else:
188            if not g.getShape()==(self.domain.getDim(),):
189               raise ValueError,"illegal shape of g"
190            self.__g=g
191          if location_of_fixed_pressure!=None:
192               self.__pde_p.setValue(q=location_of_fixed_pressure)
193               #self.__pde_l.setValue(q=location_of_fixed_pressure)
194          if location_of_fixed_flux!=None:
195               self.__pde_k.setValue(q=location_of_fixed_flux)
196                
197          if permeability!=None:
198         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
199         if perm.getRank()==0:
200            perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
201            perm=perm*util.kronecker(self.domain.getDim())
202         elif perm.getRank()==2:
203            perm_inv=util.inverse(perm)
204         else:
205            raise ValueError,"illegal rank of permeability."
206    
207         self.__permeability=perm
208         self.__permeability_inv=perm_inv
209         self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
210         #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
211         if  self.solveForFlux:
212            self.__pde_k.setValue(D=self.__permeability_inv)
213         else:
214            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
215         self.__pde_p.setValue(A=self.__permeability)
216         #self.__pde_l.setValue(D=1/self.__l)
217             #self.__pde_l.setValue(A=self.__permeability)
218    
219       def __applWeight(self, v, f=None):
220          # solves L p = f-Dv with p = 0
221          if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
222          if f == None:
223         return -util.div(v)*self.__l
224          else:
225         return (f-util.div(v))*self.__l
226          # if f == None:
227          #      self.__pde_l.setValue(Y=-util.div(v))  
228          # else:
229          #      return (f-util.div(v))/self.__l
230          # return self.__pde_l.getSolution()
231          
232       def __getPressure(self, v, p0, g=None):
233          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
234          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
235          if g == None:
236         self.__pde_p.setValue(X=-v, r=p0)
237          else:
238         self.__pde_p.setValue(X=g-v, r=p0)      
239          p=self.__pde_p.getSolution()
240          return p
241    
242       def __Aprod_v(self,dv):
243          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
244          dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
246          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247          return ArithmeticTuple(a,b,-dp)
248    
249       def __Msolve_PCG_v(self,r):
250          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())
253          # self.__pde_p.getOperator().saveMM("prec.mm")
254          return self.__pde_k.getSolution()
255    
256       def __inner_PCG_v(self,v,r):
257          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
258          
259       def __Aprod_p(self,dp):
260          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261          Gdp=util.grad(dp)
262          self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())
263          du=self.__pde_k.getSolution()
264          # self.__pde_v.getOperator().saveMM("proj.mm")
265          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
266    
267       def __getFlux(self,p, v0, f=None, g=None):
268          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
269          if f!=None:
270         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
271          self.__pde_k.setValue(r=v0)
272          g2=util.tensor_mult(self.__permeability_inv,g)
273          if p == None:
274         self.__pde_k.setValue(Y=g2)
275          else:
276         self.__pde_k.setValue(Y=g2-util.grad(p))
277          return self.__pde_k.getSolution()  
278          
279          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
280       def __Msolve_PCG_p(self,r):
281          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282          self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())
283          # self.__pde_p.getOperator().saveMM("prec.mm")
284          return self.__pde_p.getSolution()
285            
286       def __inner_PCG_p(self,p,r):
287           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288    
289       def __L2(self,v):
290          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291    
292       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293          """
294          solves the problem.
295          
296          The iteration is terminated if the residual norm is less then self.getTolerance().
297    
298          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
299          :type u0: vector value on the domain (e.g. `Data`).
300          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
301          :type p0: scalar value on the domain (e.g. `Data`).
302          :param verbose: if set some information on iteration progress are printed
303          :type verbose: ``bool``
304          :return: flux and pressure
305          :rtype: ``tuple`` of `Data`.
306    
307          :note: The problem is solved as a least squares form
308          *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
309          *G^*u+G^* K Gp=G^*g*
310    
311          where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
312          """
313          self.verbose=verbose
314          rtol=self.getTolerance()
315          atol=self.getAbsoluteTolerance()
316          self.setSubProblemTolerance()
317          num_corrections=0
318          converged=False
319          norm_r=None
320          
321          # Eliminate the hydrostatic pressure:
322          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
323          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
324          p0=self.__pde_p.getSolution()
325          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
326          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
327    
328          p=p0*0
329          if self.solveForFlux:
330         v=u0.copy()
331          else:
332         v=self.__getFlux(p, u0, f=self.__f, g=g2)
333    
334          while not converged and norm_g2 > 0:
335         Gp=util.grad(p)
336         KGp=util.tensor_mult(self.__permeability,Gp)
337         if self.verbose:
338            def_p=g2-(v+KGp)
339            def_v=self.__f-util.div(v)
340            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
341            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
342            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
343            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
344            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
345         ATOL=atol+rtol*norm_g2
346         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
347         if norm_r == None or norm_r>ATOL:
348            if num_corrections>max_num_corrections:
349               raise ValueError,"maximum number of correction steps reached."
350          
351            if self.solveForFlux:
352               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
353               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
354                   self.__Aprod_v,
355                   v,
356                   self.__Msolve_PCG_v,
357                   self.__inner_PCG_v,
358                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
359               p=r[2]
360            else:
361               # initial residual is r=G^*(g2-KGp - v)
362               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
363                     self.__Aprod_p,
364                     p,
365                     self.__Msolve_PCG_p,
366                     self.__inner_PCG_p,
367                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
368               v=r[1]
369            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
370            num_corrections+=1
371         else:
372            if self.verbose: print "DarcyFlux: stopping criterium reached."
373            converged=True
374          return v,p+p0
375       def setTolerance(self,rtol=1e-4):
376          """
377          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
378    
379          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
380          
381          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
382          
383          :param rtol: relative tolerance for the pressure
384          :type rtol: non-negative ``float``
385          """
386          if rtol<0:
387         raise ValueError,"Relative tolerance needs to be non-negative."
388          self.__rtol=rtol
389       def getTolerance(self):
390          """
391          returns the relative tolerance
392          :return: current relative tolerance
393          :rtype: ``float``
394          """
395          return self.__rtol
396    
397       def setAbsoluteTolerance(self,atol=0.):
398          """
399          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
400          
401          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
402    
403    
404          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
405    
406          :param atol: absolute tolerance for the pressure
407          :type atol: non-negative ``float``
408          """
409          if atol<0:
410         raise ValueError,"Absolute tolerance needs to be non-negative."
411          self.__atol=atol
412       def getAbsoluteTolerance(self):
413          """
414          returns the absolute tolerance
415          :return: current absolute tolerance
416          :rtype: ``float``
417          """
418          return self.__atol
419       def getSubProblemTolerance(self):
420          """
421          Returns a suitable subtolerance
422          :type: ``float``
423          """
424          return max(util.EPSILON**(0.75),self.getTolerance()**2)
425    
426       def setSubProblemTolerance(self):
427          """
428          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
429          """
430          if self.__adaptSubTolerance:
431         sub_tol=self.getSubProblemTolerance()
432         self.getSolverOptionsFlux().setTolerance(sub_tol)
433         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
434         self.getSolverOptionsPressure().setTolerance(sub_tol)
435         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
436         self.getSolverOptionsWeighting().setTolerance(sub_tol)
437         self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
438         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
439    
440    
441    class DarcyFlowOld(object):
442        """
443        solves the problem
444    
445        *u_i+k_{ij}*p_{,j} = g_i*
446        *u_{i,i} = f*
447    
448        where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
449    
450        :note: The problem is solved in a least squares formulation.
451        """
452    
453        def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
454            """
455            initializes the Darcy flux problem
456            :param domain: domain of the problem
457            :type domain: `Domain`
458        :param useReduced: uses reduced oreder on flux and pressure
459        :type useReduced: ``bool``
460        :param adaptSubTolerance: switches on automatic subtolerance selection
461        :type adaptSubTolerance: ``bool``  
462            """
463            self.domain=domain
464            if weight == None:
465               s=self.domain.getSize()
466               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
467               # self.__l=(3.*util.longestEdge(self.domain))**2
468               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
469            else:
470               self.__l=weight
471            self.__pde_v=LinearPDESystem(domain)
472            if useReduced: self.__pde_v.setReducedOrderOn()
473            self.__pde_v.setSymmetryOn()
474            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
475            self.__pde_p=LinearSinglePDE(domain)
476            self.__pde_p.setSymmetryOn()
477            if useReduced: self.__pde_p.setReducedOrderOn()
478            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
479            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
480            self.setTolerance()
481            self.setAbsoluteTolerance()
482        self.__adaptSubTolerance=adaptSubTolerance
483        self.verbose=False
484        def getSolverOptionsFlux(self):
485        """
486        Returns the solver options used to solve the flux problems
487        
488        *(I+D^*D)u=F*
489        
490        :return: `SolverOptions`
491        """
492        return self.__pde_v.getSolverOptions()
493        def setSolverOptionsFlux(self, options=None):
494        """
495        Sets the solver options used to solve the flux problems
496        
497        *(I+D^*D)u=F*
498        
499        If ``options`` is not present, the options are reset to default
500        :param options: `SolverOptions`
501        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
502        """
503        return self.__pde_v.setSolverOptions(options)
504        def getSolverOptionsPressure(self):
505        """
506        Returns the solver options used to solve the pressure problems
507        
508        *(Q^*Q)p=Q^*G*
509        
510        :return: `SolverOptions`
511        """
512        return self.__pde_p.getSolverOptions()
513        def setSolverOptionsPressure(self, options=None):
514        """
515        Sets the solver options used to solve the pressure problems
516        
517        *(Q^*Q)p=Q^*G*
518        
519        If ``options`` is not present, the options are reset to default
520        :param options: `SolverOptions`
521        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
522        """
523        return self.__pde_p.setSolverOptions(options)
524    
525        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
526            """
527            assigns values to model parameters
528    
529            :param f: volumetic sources/sinks
530            :type f: scalar value on the domain (e.g. `Data`)
531            :param g: flux sources/sinks
532            :type g: vector values on the domain (e.g. `Data`)
533            :param location_of_fixed_pressure: mask for locations where pressure is fixed
534            :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
535            :param location_of_fixed_flux:  mask for locations where flux is fixed.
536            :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
537            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
538                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
539                                 ``v`` on the main diagonal is used.
540            :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
541    
542            :note: the values of parameters which are not set by calling ``setValue`` are not altered.
543            :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
544                   or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
545                   is along the *x_i* axis.
546            """
547            if f !=None:
548               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
549               if f.isEmpty():
550                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
551               else:
552                   if f.getRank()>0: raise ValueError,"illegal rank of f."
553               self.__f=f
554            if g !=None:
555               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
556               if g.isEmpty():
557                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
558               else:
559                 if not g.getShape()==(self.domain.getDim(),):
560                   raise ValueError,"illegal shape of g"
561               self.__g=g
562    
563            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
564            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
565    
566            if permeability!=None:
567               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
568               if perm.getRank()==0:
569                   perm=perm*util.kronecker(self.domain.getDim())
570               elif perm.getRank()==1:
571                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
572                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
573               elif perm.getRank()==2:
574                  pass
575               else:
576                  raise ValueError,"illegal rank of permeability."
577               self.__permeability=perm
578               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
579    
580        def setTolerance(self,rtol=1e-4):
581            """
582            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
583    
584            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
585    
586            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
587    
588            :param rtol: relative tolerance for the pressure
589            :type rtol: non-negative ``float``
590            """
591            if rtol<0:
592                raise ValueError,"Relative tolerance needs to be non-negative."
593            self.__rtol=rtol
594        def getTolerance(self):
595            """
596            returns the relative tolerance
597    
598            :return: current relative tolerance
599            :rtype: ``float``
600            """
601            return self.__rtol
602    
603        def setAbsoluteTolerance(self,atol=0.):
604            """
605            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
606    
607            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
608    
609            where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
610    
611            :param atol: absolute tolerance for the pressure
612            :type atol: non-negative ``float``
613            """
614            if atol<0:
615                raise ValueError,"Absolute tolerance needs to be non-negative."
616            self.__atol=atol
617        def getAbsoluteTolerance(self):
618           """
619           returns the absolute tolerance
620          
621           :return: current absolute tolerance
622           :rtype: ``float``
623           """
624           return self.__atol
625        def getSubProblemTolerance(self):
626        """
627        Returns a suitable subtolerance
628        @type: ``float``
629        """
630        return max(util.EPSILON**(0.75),self.getTolerance()**2)
631        def setSubProblemTolerance(self):
632             """
633             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
634             """
635         if self.__adaptSubTolerance:
636             sub_tol=self.getSubProblemTolerance()
637                 self.getSolverOptionsFlux().setTolerance(sub_tol)
638             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
639             self.getSolverOptionsPressure().setTolerance(sub_tol)
640             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
641             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
642    
643            -(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):
644             """
645             solves the problem.
646    
647             The iteration is terminated if the residual norm is less then self.getTolerance().
648    
649             :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
650             :type u0: vector value on the domain (e.g. `Data`).
651             :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
652             :type p0: scalar value on the domain (e.g. `Data`).
653             :param verbose: if set some information on iteration progress are printed
654             :type verbose: ``bool``
655             :return: flux and pressure
656             :rtype: ``tuple`` of `Data`.
657    
658             :note: The problem is solved as a least squares form
659    
660             *(I+D^*D)u+Qp=D^*f+g*
661             *Q^*u+Q^*Qp=Q^*g*
662    
663             where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
664             We eliminate the flux form the problem by setting
665    
666             *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
667    
668             form the first equation. Inserted into the second equation we get
669    
670             *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
671    
672             which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
673             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
674             """
675             self.verbose=verbose
676             rtol=self.getTolerance()
677             atol=self.getAbsoluteTolerance()
678         self.setSubProblemTolerance()
679             num_corrections=0
680             converged=False
681             p=p0
682             norm_r=None
683             while not converged:
684                   v=self.getFlux(p, fixed_flux=u0)
685                   Qp=self.__Q(p)
686                   norm_v=self.__L2(v)
687                   norm_Qp=self.__L2(Qp)
688                   if norm_v == 0.:
689                      if norm_Qp == 0.:
690                         return v,p
691                      else:
692                        fac=norm_Qp
693                   else:
694                      if norm_Qp == 0.:
695                        fac=norm_v
696                      else:
697                        fac=2./(1./norm_v+1./norm_Qp)
698                   ATOL=(atol+rtol*fac)
699                   if self.verbose:
700                        print "DarcyFlux: L2 norm of v = %e."%norm_v
701                        print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
702                        print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
703                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
704                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
705                   if norm_r == None or norm_r>ATOL:
706                       if num_corrections>max_num_corrections:
707                             raise ValueError,"maximum number of correction steps reached."
708                       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)
709                       num_corrections+=1
710                   else:
711                       converged=True
712             return v,p
713        def __L2(self,v):
714             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
715    
716        def __Q(self,p):
717              return util.tensor_mult(self.__permeability,util.grad(p))
718    
719        def __Aprod(self,dp):
720              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
721              Qdp=self.__Q(dp)
722              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
723              du=self.__pde_v.getSolution()
724              # self.__pde_v.getOperator().saveMM("proj.mm")
725              return Qdp+du
726        def __inner_GMRES(self,r,s):
727             return util.integrate(util.inner(r,s))
728    
729        def __inner_PCG(self,p,r):
730             return util.integrate(util.inner(self.__Q(p), r))
731    
732        def __Msolve_PCG(self,r):
733          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
734              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
735              # self.__pde_p.getOperator().saveMM("prec.mm")
736              return self.__pde_p.getSolution()
737    
738        def getFlux(self,p=None, fixed_flux=Data()):
739            """
740            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
741            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
742            Note that ``g`` and ``f`` are used, see `setValue`.
743    
744            :param p: pressure.
745            :type p: scalar value on the domain (e.g. `Data`).
746            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
747            :type fixed_flux: vector values on the domain (e.g. `Data`).
748            :return: flux
749            :rtype: `Data`
750            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
751                   for the permeability *k_{ij}*
752            """
753        self.setSubProblemTolerance()
754            g=self.__g
755            f=self.__f
756            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
757            if p == None:
758               self.__pde_v.setValue(Y=g)
759            else:
760               self.__pde_v.setValue(Y=g-self.__Q(p))
761            return self.__pde_v.getSolution()
762    
763    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
764         """
765         solves
766    
767              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
768                  u_{i,i}=0                  u_{i,i}=0
769    
770            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
771            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
772    
773        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
774    
775        typical usage:       typical usage:
776    
777              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
778              sp.setTolerance()              sp.setTolerance()
779              sp.initialize(...)              sp.initialize(...)
780              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
781        """       """
782        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
783             """
784             initialize the Stokes Problem
785    
786             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
787             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
788             with macro elements for the pressure.
789    
790             :param domain: domain of the problem.
791             :type domain: `Domain`
792             """
793           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
794           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
795           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())
796           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
797           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)      
               
798           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
799           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
800           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
801    
802           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
803           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
804         self.__pde_proj.setValue(D=1)
805           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
806    
807        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
808          self.eta=eta           """
809          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
810      self.__pde_u.setValue(A=Data())      
811          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
812          for j in range(self.domain.getDim()):       """
813              A[i,j,j,i] += 1.       return self.__pde_u.getSolverOptions()
814              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
815      self.__pde_prec.setValue(D=1/self.eta)           """
816          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)       set the solver options for solving the equation for velocity.
817        
818        def B(self,arg):       :param options: new solver  options
819           d=util.div(arg)       :type options: `SolverOptions`
820           self.__pde_proj.setValue(Y=d)       """
821           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setSolverOptions(options)
822           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
823             """
824         returns the solver options used  solve the equation for pressure.
825         :rtype: `SolverOptions`
826         """
827         return self.__pde_prec.getSolverOptions()
828         def setSolverOptionsPressure(self, options=None):
829             """
830         set the solver options for solving the equation for pressure.
831         :param options: new solver  options
832         :type options: `SolverOptions`
833         """
834         self.__pde_prec.setSolverOptions(options)
835    
836         def setSolverOptionsDiv(self, options=None):
837             """
838         set the solver options for solving the equation to project the divergence of
839         the velocity onto the function space of presure.
840        
841         :param options: new solver options
842         :type options: `SolverOptions`
843         """
844         self.__pde_proj.setSolverOptions(options)
845         def getSolverOptionsDiv(self):
846             """
847         returns the solver options for solving the equation to project the divergence of
848         the velocity onto the function space of presure.
849        
850         :rtype: `SolverOptions`
851         """
852         return self.__pde_proj.getSolverOptions()
853    
854         def updateStokesEquation(self, v, p):
855             """
856             updates the Stokes equation to consider dependencies from ``v`` and ``p``
857             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
858             """
859             pass
860         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
861            """
862            assigns new values to the model parameters.
863    
864            :param f: external force
865            :type f: `Vector` object in `FunctionSpace` `Function` or similar
866            :param fixed_u_mask: mask of locations with fixed velocity.
867            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
868            :param eta: viscosity
869            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
870            :param surface_stress: normal surface stress
871            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
872            :param stress: initial stress
873        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
874            """
875            if eta !=None:
876                k=util.kronecker(self.domain.getDim())
877                kk=util.outer(k,k)
878                self.eta=util.interpolate(eta, Function(self.domain))
879            self.__pde_prec.setValue(D=1/self.eta)
880                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
881            if restoration_factor!=None:
882                n=self.domain.getNormal()
883                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
884            if fixed_u_mask!=None:
885                self.__pde_u.setValue(q=fixed_u_mask)
886            if f!=None: self.__f=f
887            if surface_stress!=None: self.__surface_stress=surface_stress
888            if stress!=None: self.__stress=stress
889    
890         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
891            """
892            assigns values to the model parameters
893    
894            :param f: external force
895            :type f: `Vector` object in `FunctionSpace` `Function` or similar
896            :param fixed_u_mask: mask of locations with fixed velocity.
897            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
898            :param eta: viscosity
899            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
900            :param surface_stress: normal surface stress
901            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
902            :param stress: initial stress
903        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
904            """
905            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
906    
907         def Bv(self,v,tol):
908             """
909             returns inner product of element p and div(v)
910    
911             :param v: a residual
912             :return: inner product of element p and div(v)
913             :rtype: ``float``
914             """
915             self.__pde_proj.setValue(Y=-util.div(v))
916         self.getSolverOptionsDiv().setTolerance(tol)
917         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
918             out=self.__pde_proj.getSolution()
919             return out
920    
921         def inner_pBv(self,p,Bv):
922             """
923             returns inner product of element p and Bv=-div(v)
924    
925             :param p: a pressure increment
926             :param Bv: a residual
927             :return: inner product of element p and Bv=-div(v)
928             :rtype: ``float``
929             """
930             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
931    
932        def inner(self,p0,p1):       def inner_p(self,p0,p1):
933             """
934             Returns inner product of p0 and p1
935    
936             :param p0: a pressure
937             :param p1: a pressure
938             :return: inner product of p0 and p1
939             :rtype: ``float``
940             """
941           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
942           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
943           return util.integrate(s0*s1)           return util.integrate(s0*s1)
944    
945        def inner_a(self,a0,a1):       def norm_v(self,v):
946           p0=util.interpolate(a0[1],Function(self.domain))           """
947           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
          alfa=(1/self.vol)*util.integrate(p0)  
          beta=(1/self.vol)*util.integrate(p1)  
      v0=util.grad(a0[0])  
      v1=util.grad(a1[0])  
      #print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))  
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
   
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
   
       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  
           print abs(n_v_old-self.__n_v), n_v, self.getTolerance()  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():  
               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  
948    
949             :param v: a velovity
950             :return: norm of v
951             :rtype: non-negative ``float``
952             """
953             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
954    
955    
956         def getDV(self, p, v, tol):
957             """
958             return the value for v for a given p (overwrite)
959    
960             :param p: a pressure
961             :param v: a initial guess for the value v to return.
962             :return: dv given as *Adv=(f-Av-B^*p)*
963             """
964             self.updateStokesEquation(v,p)
965             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
966         self.getSolverOptionsVelocity().setTolerance(tol)
967         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
968             if self.__stress.isEmpty():
969                self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
970             else:
971                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
972             out=self.__pde_u.getSolution()
973             return  out
974    
975         def norm_Bv(self,Bv):
976            """
977            Returns Bv (overwrite).
978    
979            :rtype: equal to the type of p
980            :note: boundary conditions on p should be zero!
981            """
982            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
983    
984         def solve_AinvBt(self,p, tol):
985             """
986             Solves *Av=B^*p* with accuracy `tol`
987    
988             :param p: a pressure increment
989             :return: the solution of *Av=B^*p*
990             :note: boundary conditions on v should be zero!
991             """
992             self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993             out=self.__pde_u.getSolution()
994             return  out
995    
996         def solve_prec(self,Bv, tol):
997             """
998             applies preconditioner for for *BA^{-1}B^** to *Bv*
999             with accuracy `self.getSubProblemTolerance()`
1000    
1001             :param Bv: velocity increment
1002             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1003             :note: boundary conditions on p are zero.
1004             """
1005             self.__pde_prec.setValue(Y=Bv)
1006         self.getSolverOptionsPressure().setTolerance(tol)
1007         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008             out=self.__pde_prec.getSolution()
1009             return out

Legend:
Removed from v.1562  
changed lines
  Added in v.3071

  ViewVC Help
Powered by ViewVC 1.1.26