/[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 2344 by jfenwick, Mon Mar 30 02:13:58 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-2008 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-2008 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"
34    
35  from escript import *  from escript import *
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
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,useReduced=False):      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
459        :type useReduced: ``bool``
460        :param adaptSubTolerance: switches on automatic subtolerance selection
461        :type adaptSubTolerance: ``bool``  
462          """          """
463          self.domain=domain          self.domain=domain
464          self.__l=util.longestEdge(self.domain)**2          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)          self.__pde_v=LinearPDESystem(domain)
472          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
473          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
# Line 67  class DarcyFlow(object): Line 479  class DarcyFlow(object):
479          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
480          self.setTolerance()          self.setTolerance()
481          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
482          self.setSubProblemTolerance()      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):      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          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 126  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 142  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 165  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
   
     def setSubProblemTolerance(self,rtol=None):  
          """  
          Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present  
          C{self.getTolerance()**2} is used.  
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
          """  
          if rtol == None:  
               if self.getTolerance()<=0.:  
                   raise ValueError,"A positive relative tolerance must be set."  
               self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)  
          else:  
              if rtol<=0:  
                  raise ValueError,"sub-problem tolerance must be positive."  
              self.__sub_tol=max(util.EPSILON**(0.75),rtol)  
   
625      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
626           """      """
627           Returns the subproblem reduction factor.      Returns a suitable subtolerance
628        @type: ``float``
629           @return: subproblem reduction factor      """
630           @rtype: C{float}      return max(util.EPSILON**(0.75),self.getTolerance()**2)
631           """      def setSubProblemTolerance(self):
632           return self.__sub_tol           """
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      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
644           """           """
645           solves the problem.           solves the problem.
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           @param show_details:  if set information on the subiteration process are printed.           :return: flux and pressure
656           @type show_details: C{bool}           :rtype: ``tuple`` of `Data`.
          @return: flux and pressure  
          @rtype: C{tuple} of L{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 or True           self.verbose=verbose
          self.show_details= show_details and self.verbose  
676           rtol=self.getTolerance()           rtol=self.getTolerance()
677           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
678           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()       self.setSubProblemTolerance()
   
679           num_corrections=0           num_corrections=0
680           converged=False           converged=False
681           p=p0           p=p0
682           norm_r=None           norm_r=None
683           while not converged:           while not converged:
684                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)                 v=self.getFlux(p, fixed_flux=u0)
685                 Qp=self.__Q(p)                 Qp=self.__Q(p)
686                 norm_v=self.__L2(v)                 norm_v=self.__L2(v)
687                 norm_Qp=self.__L2(Qp)                 norm_Qp=self.__L2(Qp)
# Line 258  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),)
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                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
705                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
706                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
707                           raise ValueError,"maximum number of correction steps reached."                           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.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     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                     num_corrections+=1
710                 else:                 else:
711                     converged=True                     converged=True
712           return v,p           return v,p
 #  
 #                
 #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp  
 #               #===========================================================================  
 #               norm_r_hat=self.__L2(r_hat)  
 #               norm_v=self.__L2(v)  
 #               norm_g=self.__L2(g)  
 #               norm_gv=self.__L2(g-v)  
 #               norm_Qp=self.__L2(Qp)  
 #               norm_gQp=self.__L2(g-Qp)  
 #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))  
 #               fac=min(norm_v,norm_Qp,norm_gv)  
 #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat  
 #               if r!=None:  
 #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat  
 #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()  
 #                   self.setSubProblemTolerance(sub_tol)  
 #                   print "subtol_new=",self.getSubProblemTolerance()  
 #               print "norm_v = ",norm_v  
 #               print "norm_gv = ",norm_gv  
 #               print "norm_Qp = ",norm_Qp  
 #               print "norm_gQp = ",norm_gQp  
 #               print "norm_g = ",norm_g  
 #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)  
 #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)  
 #               if fac == 0:  
 #                   if self.verbose: print "DarcyFlux: trivial case!"  
 #                   return v,p  
 #               #===============================================================================  
 #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))  
 #               # norm_Qp=self.__L2(Qp)  
 #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp  
 #  
 #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))  
 #               if self.verbose:  
 #                   print "DarcyFlux: residual = %e"%norm_r_hat  
 #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
 #               if norm_r_hat <= ATOL:  
 #                   print "DarcyFlux: iteration finalized."  
 #                   converged=True  
 #               else:  
 #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)  
 #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #               print "norm_r =",norm_r  
 #         return v,p  
713      def __L2(self,v):      def __L2(self,v):
714           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
715    
# Line 323  class DarcyFlow(object): Line 717  class DarcyFlow(object):
717            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
718    
719      def __Aprod(self,dp):      def __Aprod(self,dp):
720            self.__pde_v.setTolerance(self.getSubProblemTolerance())            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
           if self.show_details: print "DarcyFlux: Applying operator"  
721            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
722            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
723            du=self.__pde_v.getSolution(verbose=self.show_details)            du=self.__pde_v.getSolution()
724              # self.__pde_v.getOperator().saveMM("proj.mm")
725            return Qdp+du            return Qdp+du
726      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
727           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 336  class DarcyFlow(object): Line 730  class DarcyFlow(object):
730           return util.integrate(util.inner(self.__Q(p), r))           return util.integrate(util.inner(self.__Q(p), r))
731    
732      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
733            self.__pde_p.setTolerance(self.getSubProblemTolerance())        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
           if self.show_details: print "DarcyFlux: Applying preconditioner"  
734            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
735            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
736              return self.__pde_p.getSolution()
737    
738      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):      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.__pde_v.setTolerance(self.getSubProblemTolerance())      self.setSubProblemTolerance()
754          g=self.__g          g=self.__g
755          f=self.__f          f=self.__f
756          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
# Line 366  class DarcyFlow(object): Line 758  class DarcyFlow(object):
758             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
759          else:          else:
760             self.__pde_v.setValue(Y=g-self.__Q(p))             self.__pde_v.setValue(Y=g-self.__Q(p))
761          return self.__pde_v.getSolution(verbose=show_details)          return self.__pde_v.getSolution()
762    
763  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
764       """       """
# Line 391  class StokesProblemCartesian(Homogeneous Line 783  class StokesProblemCartesian(Homogeneous
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           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           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(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
798           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
799           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
800           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
801    
802       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):           self.__pde_proj=LinearPDE(domain)
803             self.__pde_proj.setReducedOrderOn()
804         self.__pde_proj.setValue(D=1)
805             self.__pde_proj.setSymmetryOn()
806    
807         def getSolverOptionsVelocity(self):
808             """
809         returns the solver options used  solve the equation for velocity.
810        
811         :rtype: `SolverOptions`
812         """
813         return self.__pde_u.getSolverOptions()
814         def setSolverOptionsVelocity(self, options=None):
815             """
816         set the solver options for solving the equation for velocity.
817        
818         :param options: new solver  options
819         :type options: `SolverOptions`
820         """
821             self.__pde_u.setSolverOptions(options)
822         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          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.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
906    
907       def inner_pBv(self,p,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.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(-p*util.div(v))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
931    
932       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
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.setTolerance(self.getSubProblemTolerance())           self.updateStokesEquation(v,p)
965           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           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(verbose=self.show_details)           out=self.__pde_u.getSolution()
973           return  out           return  out
974    
975         def norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
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.div(v)**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.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           out=self.__pde_u.getSolution()
          out=self.__pde_u.getSolution(verbose=self.show_details)  
994           return  out           return  out
995    
996       def solve_precB(self,v):       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()} (overwrite).           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=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
1006           self.__pde_prec.setTolerance(self.getSubProblemTolerance())       self.getSolverOptionsPressure().setTolerance(tol)
1007           return self.__pde_prec.getSolution(verbose=self.show_details)       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008             out=self.__pde_prec.getSolution()
1009             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26