/[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 2620 by gross, Thu Aug 20 06:24:00 2009 UTC revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
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"
# Line 37  from linearPDEs import LinearPDE, Linear Line 38  from linearPDEs import LinearPDE, Linear
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  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          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.5),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      solves the problem
444    
445      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
446      M{u_{i,i} = f}      *u_{i,i} = f*
447    
448      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      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.      :note: The problem is solved in a least squares formulation.
451      """      """
452    
453      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
454          """          """
455          initializes the Darcy flux problem          initializes the Darcy flux problem
456          @param domain: domain of the problem          :param domain: domain of the problem
457          @type domain: L{Domain}          :type domain: `Domain`
458      @param useReduced: uses reduced oreder on flux and pressure      :param useReduced: uses reduced oreder on flux and pressure
459      @type useReduced: C{bool}      :type useReduced: ``bool``
460      @param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
461      @type adaptSubTolerance: C{bool}          :type adaptSubTolerance: ``bool``  
462          """          """
463          self.domain=domain          self.domain=domain
464          if weight == None:          if weight == None:
465             s=self.domain.getSize()             s=self.domain.getSize()
466             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
467             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l=(3.*util.longestEdge(self.domain))**2
468             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
469          else:          else:
470             self.__l=weight             self.__l=weight
471          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
# Line 83  class DarcyFlow(object): Line 485  class DarcyFlow(object):
485      """      """
486      Returns the solver options used to solve the flux problems      Returns the solver options used to solve the flux problems
487            
488      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
489            
490      @return: L{SolverOptions}      :return: `SolverOptions`
491      """      """
492      return self.__pde_v.getSolverOptions()      return self.__pde_v.getSolverOptions()
493      def setSolverOptionsFlux(self, options=None):      def setSolverOptionsFlux(self, options=None):
494      """      """
495      Sets the solver options used to solve the flux problems      Sets the solver options used to solve the flux problems
496            
497      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
498            
499      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
500      @param options: L{SolverOptions}      :param options: `SolverOptions`
501      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :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)      return self.__pde_v.setSolverOptions(options)
504      def getSolverOptionsPressure(self):      def getSolverOptionsPressure(self):
505      """      """
506      Returns the solver options used to solve the pressure problems      Returns the solver options used to solve the pressure problems
507            
508      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
509            
510      @return: L{SolverOptions}      :return: `SolverOptions`
511      """      """
512      return self.__pde_p.getSolverOptions()      return self.__pde_p.getSolverOptions()
513      def setSolverOptionsPressure(self, options=None):      def setSolverOptionsPressure(self, options=None):
514      """      """
515      Sets the solver options used to solve the pressure problems      Sets the solver options used to solve the pressure problems
516            
517      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
518            
519      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
520      @param options: L{SolverOptions}      :param options: `SolverOptions`
521      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :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)      return self.__pde_p.setSolverOptions(options)
524    
# Line 124  class DarcyFlow(object): Line 526  class DarcyFlow(object):
526          """          """
527          assigns values to model parameters          assigns values to model parameters
528    
529          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
530          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `Data`)
531          @param g: flux sources/sinks          :param g: flux sources/sinks
532          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `Data`)
533          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :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. L{Data})          :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.          :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. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
537          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
538                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
539                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
540          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :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 C{setValue} are not altered.          :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 (C{location_of_fixed_pressure} >0)          :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 (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
545                 is along the M{x_i} axis.                 is along the *x_i* axis.
546          """          """
547          if f !=None:          if f !=None:
548             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
# Line 177  class DarcyFlow(object): Line 579  class DarcyFlow(object):
579    
580      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
581          """          """
582          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
583    
584          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
585    
586          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}}.          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          :param rtol: relative tolerance for the pressure
589          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
590          """          """
591          if rtol<0:          if rtol<0:
592              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 193  class DarcyFlow(object): Line 595  class DarcyFlow(object):
595          """          """
596          returns the relative tolerance          returns the relative tolerance
597    
598          @return: current relative tolerance          :return: current relative tolerance
599          @rtype: C{float}          :rtype: ``float``
600          """          """
601          return self.__rtol          return self.__rtol
602    
603      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
604          """          """
605          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
606    
607          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
608    
609          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}}.          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          :param atol: absolute tolerance for the pressure
612          @type atol: non-negative C{float}          :type atol: non-negative ``float``
613          """          """
614          if atol<0:          if atol<0:
615              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 216  class DarcyFlow(object): Line 618  class DarcyFlow(object):
618         """         """
619         returns the absolute tolerance         returns the absolute tolerance
620                
621         @return: current absolute tolerance         :return: current absolute tolerance
622         @rtype: C{float}         :rtype: ``float``
623         """         """
624         return self.__atol         return self.__atol
625      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
626      """      """
627      Returns a suitable subtolerance      Returns a suitable subtolerance
628      @type: C{float}      @type: ``float``
629      """      """
630      return max(util.EPSILON**(0.75),self.getTolerance()**2)      return max(util.EPSILON**(0.75),self.getTolerance()**2)
631      def setSubProblemTolerance(self):      def setSubProblemTolerance(self):
# Line 244  class DarcyFlow(object): Line 646  class DarcyFlow(object):
646    
647           The iteration is terminated if the residual norm is less then self.getTolerance().           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 C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           :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. L{Data}).           :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 C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           :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. L{Data}).           :type p0: scalar value on the domain (e.g. `Data`).
653           @param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
654           @type verbose: C{bool}           :type verbose: ``bool``
655           @return: flux and pressure           :return: flux and pressure
656           @rtype: C{tuple} of L{Data}.           :rtype: ``tuple`` of `Data`.
657    
658           @note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
659    
660           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
661           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
662    
663           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           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           We eliminate the flux form the problem by setting
665    
666           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *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           form the first equation. Inserted into the second equation we get
669    
670           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           *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 M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
673           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
674           """           """
675           self.verbose=verbose           self.verbose=verbose
676           rtol=self.getTolerance()           rtol=self.getTolerance()
677           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
678       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
679           num_corrections=0           num_corrections=0
680           converged=False           converged=False
681           p=p0           p=p0
# Line 297  class DarcyFlow(object): Line 698  class DarcyFlow(object):
698                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
699                 if self.verbose:                 if self.verbose:
700                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
701                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      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),)                      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),)                      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                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
# Line 336  class DarcyFlow(object): Line 737  class DarcyFlow(object):
737    
738      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=Data()):
739          """          """
740          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
741          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
742          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
743    
744          @param p: pressure.          :param p: pressure.
745          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `Data`).
746          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          :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. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `Data`).
748          @param tol: relative tolerance to be used.          :return: flux
749          @type tol: positive C{float}.          :rtype: `Data`
750          @return: flux          :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          @rtype: L{Data}                 for the permeability *k_{ij}*
         @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}}  
                for the permeability M{k_{ij}}  
752          """          """
753      self.setSubProblemTolerance()      self.setSubProblemTolerance()
754          g=self.__g          g=self.__g
# Line 380  class StokesProblemCartesian(Homogeneous Line 779  class StokesProblemCartesian(Homogeneous
779              sp.initialize(...)              sp.initialize(...)
780              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
781       """       """
782       def __init__(self,domain,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
783           """           """
784           initialize the Stokes Problem           initialize the Stokes Problem
785    
786           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
787           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
788       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
789       @type adaptSubTolerance: C{bool}  
790           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
791             :type domain: `Domain`
792           """           """
793           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**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            
# Line 409  class StokesProblemCartesian(Homogeneous Line 808  class StokesProblemCartesian(Homogeneous
808           """           """
809       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
810            
811       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
812       """       """
813       return self.__pde_u.getSolverOptions()       return self.__pde_u.getSolverOptions()
814       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
815           """           """
816       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
817            
818       @param options: new solver  options       :param options: new solver  options
819       @type options: L{SolverOptions}       :type options: `SolverOptions`
820       """       """
821           self.__pde_u.setSolverOptions(options)           self.__pde_u.setSolverOptions(options)
822       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
823           """           """
824       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
825       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
826       """       """
827       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
828       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
829           """           """
830       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
831       @param options: new solver  options       :param options: new solver  options
832       @type options: L{SolverOptions}       :type options: `SolverOptions`
833       """       """
834       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
835    
# Line 439  class StokesProblemCartesian(Homogeneous Line 838  class StokesProblemCartesian(Homogeneous
838       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
839       the velocity onto the function space of presure.       the velocity onto the function space of presure.
840            
841       @param options: new solver options       :param options: new solver options
842       @type options: L{SolverOptions}       :type options: `SolverOptions`
843       """       """
844       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
845       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
846           """           """
847       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
848       the velocity onto the function space of presure.       the velocity onto the function space of presure.
849            
850       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
851       """       """
852       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
853       def setSubProblemTolerance(self):  
854         def updateStokesEquation(self, v, p):
855           """           """
856       Updates the tolerance for subproblems           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       if self.adaptSubTolerance():           pass
860               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
861           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
862           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
863           self.getSolverOptionsPressure().setTolerance(sub_tol)  
864           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
865           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
866           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :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):       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          assigns values to the model parameters
893    
894          @param f: external force          :param f: external force
895          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
896          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
897          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
898          @param eta: viscosity          :param eta: viscosity
899          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
900          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
901          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
902          @param stress: initial stress          :param stress: initial stress
903      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
904          @note: All values needs to be set.          """
905          """          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         self.eta=eta  
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
         n=self.domain.getNormal()  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
906    
907       def Bv(self,v):       def Bv(self,v,tol):
908           """           """
909           returns inner product of element p and div(v)           returns inner product of element p and div(v)
910    
911           @param p: a pressure increment           :param v: a residual
912           @param v: a residual           :return: inner product of element p and div(v)
913           @return: inner product of element p and div(v)           :rtype: ``float``
914           @rtype: C{float}           """
915           """           self.__pde_proj.setValue(Y=-util.div(v))
916           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
917           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
918             out=self.__pde_proj.getSolution()
919             return out
920    
921       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
922           """           """
923           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
924    
925           @param p: a pressure increment           :param p: a pressure increment
926           @param v: a residual           :param Bv: a residual
927           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
928           @rtype: C{float}           :rtype: ``float``
929           """           """
930           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
931    
# Line 522  class StokesProblemCartesian(Homogeneous Line 933  class StokesProblemCartesian(Homogeneous
933           """           """
934           Returns inner product of p0 and p1           Returns inner product of p0 and p1
935    
936           @param p0: a pressure           :param p0: a pressure
937           @param p1: a pressure           :param p1: a pressure
938           @return: inner product of p0 and p1           :return: inner product of p0 and p1
939           @rtype: C{float}           :rtype: ``float``
940           """           """
941           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
942           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
943           return util.integrate(s0*s1)           return util.integrate(s0*s1)
944    
945       def norm_v(self,v):       def norm_v(self,v):
946           """           """
947           returns the norm of v           returns the norm of v
948    
949           @param v: a velovity           :param v: a velovity
950           @return: norm of v           :return: norm of v
951           @rtype: non-negative C{float}           :rtype: non-negative ``float``
952           """           """
953           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
954    
955    
956       def getV(self, p, v0):       def getDV(self, p, v, tol):
957           """           """
958           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
959    
960           @param p: a pressure           :param p: a pressure
961           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
962           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
963           """           """
964           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           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():           if self.__stress.isEmpty():
969              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
970           else:           else:
971              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              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()           out=self.__pde_u.getSolution()
973           return  out           return  out
974    
# Line 561  class StokesProblemCartesian(Homogeneous Line 976  class StokesProblemCartesian(Homogeneous
976          """          """
977          Returns Bv (overwrite).          Returns Bv (overwrite).
978    
979          @rtype: equal to the type of p          :rtype: equal to the type of p
980          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
981          """          """
982          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
983    
984       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
985           """           """
986           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
987    
988           @param p: a pressure increment           :param p: a pressure increment
989           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
990           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
991           """           """
992           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
994           return  out           return  out
995    
996       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
997           """           """
998           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
999           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
1000    
1001           @param v: velocity increment           :param Bv: velocity increment
1002           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1003           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
1004           """           """
1005           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
1006           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
1007         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008             out=self.__pde_prec.getSolution()
1009             return out

Legend:
Removed from v.2620  
changed lines
  Added in v.3074

  ViewVC Help
Powered by ViewVC 1.1.26