/[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 2169 by caltinay, Wed Dec 17 03:08:58 2008 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"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
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  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      Represents and solves the problem     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      M{u_i+k_{ij}*p_{,j} = g_i}     def __L2(self,v):
290          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291    
292      M{u_{i,i} = f}     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 M{p} represents the pressure and M{u} the Darcy flux. M{k} represents        where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
312      the permeability.        """
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      @note: The problem is solved in a least squares formulation.        *|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
444    
445      def __init__(self, domain):      *u_i+k_{ij}*p_{,j} = g_i*
446          """      *u_{i,i} = f*
447          Initializes the Darcy flux problem.  
448        where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
449    
450        :note: The problem is solved in a least squares formulation.
451        """
452    
453          @param domain: domain of the problem      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
454          @type domain: L{Domain}          """
455            initializes the Darcy flux problem
456            :param domain: domain of the problem
457            :type domain: `Domain`
458        :param useReduced: uses reduced oreder on flux and pressure
459        :type useReduced: ``bool``
460        :param adaptSubTolerance: switches on automatic subtolerance selection
461        :type adaptSubTolerance: ``bool``  
462          """          """
463          self.domain=domain          self.domain=domain
464            if weight == None:
465               s=self.domain.getSize()
466               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
467               # self.__l=(3.*util.longestEdge(self.domain))**2
468               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
469            else:
470               self.__l=weight
471          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
472          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          if useReduced: self.__pde_v.setReducedOrderOn()
473          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
474            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
475          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
476          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
477            if useReduced: self.__pde_p.setReducedOrderOn()
478          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
479          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
480            self.setTolerance()
481            self.setAbsoluteTolerance()
482        self.__adaptSubTolerance=adaptSubTolerance
483        self.verbose=False
484        def getSolverOptionsFlux(self):
485        """
486        Returns the solver options used to solve the flux problems
487        
488        *(I+D^*D)u=F*
489        
490        :return: `SolverOptions`
491        """
492        return self.__pde_v.getSolverOptions()
493        def setSolverOptionsFlux(self, options=None):
494        """
495        Sets the solver options used to solve the flux problems
496        
497        *(I+D^*D)u=F*
498        
499        If ``options`` is not present, the options are reset to default
500        :param options: `SolverOptions`
501        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
502        """
503        return self.__pde_v.setSolverOptions(options)
504        def getSolverOptionsPressure(self):
505        """
506        Returns the solver options used to solve the pressure problems
507        
508        *(Q^*Q)p=Q^*G*
509        
510        :return: `SolverOptions`
511        """
512        return self.__pde_p.getSolverOptions()
513        def setSolverOptionsPressure(self, options=None):
514        """
515        Sets the solver options used to solve the pressure problems
516        
517        *(Q^*Q)p=Q^*G*
518        
519        If ``options`` is not present, the options are reset to default
520        :param options: `SolverOptions`
521        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
522        """
523        return self.__pde_p.setSolverOptions(options)
524    
525      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      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: volumetric 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 value 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 value 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          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
538                               tensor with C{s} on the main diagonal is used. If                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
539                               vector C{v} is given the tensor with C{v} on the                               ``v`` on the main diagonal is used.
540                               main diagonal is used.          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
541          @type permeability: scalar, vector or tensor values on the domain, e.g.  
542                              L{Data}          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
543            :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
544          @note: the values of parameters which are not set by calling                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
545                 C{setValue} are not altered                 is along the *x_i* axis.
         @note: at any point on the boundary of the domain the pressure  
                (C{location_of_fixed_pressure}) >0 or the normal component of  
                the flux (C{location_of_fixed_flux[i]}) >0 if the direction of  
                the normal is along the M{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 98  class DarcyFlow(object): Line 550  class DarcyFlow(object):
550                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
551             else:             else:
552                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
553             self.f=f             self.__f=f
554          if g !=None:          if g !=None:
555             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
556             if g.isEmpty():             if g.isEmpty():
# Line 125  class DarcyFlow(object): Line 577  class DarcyFlow(object):
577             self.__permeability=perm             self.__permeability=perm
578             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
579    
580        def setTolerance(self,rtol=1e-4):
581            """
582            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
583    
584            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
585    
586      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
587    
588            :param rtol: relative tolerance for the pressure
589            :type rtol: non-negative ``float``
590            """
591            if rtol<0:
592                raise ValueError,"Relative tolerance needs to be non-negative."
593            self.__rtol=rtol
594        def getTolerance(self):
595          """          """
596          Returns the flux for a given pressure C{p}.          returns the relative tolerance
597    
598          The flux is equal to C{fixed_flux} on locations where          :return: current relative tolerance
599          C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}          :rtype: ``float``
600          and C{f} are used.          """
601            return self.__rtol
         @param p: pressure  
         @type p: scalar value on the domain, e.g. L{Data}  
         @param fixed_flux: flux on the locations of the domain marked by  
                            C{location_of_fixed_flux}  
         @type fixed_flux: vector values on the domain, e.g. L{Data}  
         @param tol: relative tolerance to be used  
         @type tol: positive float  
         @return: flux  
         @rtype: L{Data}  
         @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}}  
         """  
         self.__pde_v.setTolerance(tol)  
         self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)  
         return self.__pde_v.getSolution(verbose=show_details)  
   
     def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):  
         """  
         Solves the problem.  
   
         The iteration is terminated if the error in the pressure is less than  
         M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand  
         side (see escript user's guide for details).  
   
         @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.  
         @type u0: vector value on the domain, e.g. L{Data}  
         @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.  
         @type p0: scalar value on the domain, e.g. L{Data}  
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @param sub_rtol: tolerance to be used in the sub iteration. It is  
                          recommended that M{sub_rtol<rtol*5.e-3}  
         @type sub_rtol: positive-negative C{float}  
         @param verbose: if True information on iteration progress is printed  
         @type verbose: C{bool}  
         @param show_details: if True information on the sub-iteration process  
                              is printed  
         @type show_details: C{bool}  
         @return: flux and pressure  
         @rtype: C{tuple} of L{Data}  
   
         @note: the problem is solved in a least squares formulation:  
   
         M{(I+D^*D)u+Qp=D^*f+g}  
   
         M{Q^*u+Q^*Qp=Q^*g}  
   
         where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the  
         permeability M{k_{ij}}. We eliminate the flux from the problem by  
         setting  
   
         M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}  
   
         from the first equation. Inserted into the second equation we get  
   
         M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}  
         on C{location_of_fixed_pressure}  
   
         which is solved using the PCG method (precondition is M{Q^*Q}).  
         In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}  
         need to be solved using a sub-iteration scheme.  
         """  
         self.verbose=verbose  
         self.show_details= show_details and self.verbose  
         self.__pde_v.setTolerance(sub_rtol)  
         self.__pde_p.setTolerance(sub_rtol)  
         u2=u0*self.__pde_v.getCoefficient("q")  
         #  
         # first the reference velocity is calculated from  
         #  
         #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)  
         #  
         self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)  
         u_ref=self.__pde_v.getSolution(verbose=show_details)  
         if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)  
         self.__pde_v.setValue(r=Data())  
         #  
         #   and then we calculate a reference pressure  
         #  
         #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)  
         #  
         self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)  
         p_ref=self.__pde_p.getSolution(verbose=self.show_details)  
         if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)  
         self.__pde_p.setValue(r=Data())  
         #  
         #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref  
         #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref  
         #  
         #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du  
         #  
         #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref  
         #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref  
         #  
         #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)  
         #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)  
         #  
         #  we use the (du,Qdp) to represent the resudual  
         #  Q^*Q is a preconditioner  
         #  
         #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>  
         #  
         Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))  
         norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))  
         ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)  
         if not ATOL>0:  
             raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)  
         if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)  
         #  
         #   caclulate the initial residual  
         #  
         self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())  
         du=self.__pde_v.getSolution(verbose=show_details)  
         r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)  
         dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
         util.saveVTK("d.vtu",p=dp,p_ref=p_ref)  
         return u_ref+r[1],dp  
   
     def __Aprod_PCG(self,p):  
         if self.show_details: print "DarcyFlux: Applying operator"  
         Qp=util.tensor_mult(self.__permeability,util.grad(p))  
         self.__pde_v.setValue(Y=Qp,X=Data())  
         w=self.__pde_v.getSolution(verbose=self.show_details)  
         return ArithmeticTuple(-Qp,w)  
602    
603        def setAbsoluteTolerance(self,atol=0.):
604            """
605            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
606    
607            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
608    
609            where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
610    
611            :param atol: absolute tolerance for the pressure
612            :type atol: non-negative ``float``
613            """
614            if atol<0:
615                raise ValueError,"Absolute tolerance needs to be non-negative."
616            self.__atol=atol
617        def getAbsoluteTolerance(self):
618           """
619           returns the absolute tolerance
620          
621           :return: current absolute tolerance
622           :rtype: ``float``
623           """
624           return self.__atol
625        def getSubProblemTolerance(self):
626        """
627        Returns a suitable subtolerance
628        @type: ``float``
629        """
630        return max(util.EPSILON**(0.75),self.getTolerance()**2)
631        def setSubProblemTolerance(self):
632             """
633             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
634             """
635         if self.__adaptSubTolerance:
636             sub_tol=self.getSubProblemTolerance()
637                 self.getSolverOptionsFlux().setTolerance(sub_tol)
638             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
639             self.getSolverOptionsPressure().setTolerance(sub_tol)
640             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
641             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
642    
643        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
644             """
645             solves the problem.
646    
647             The iteration is terminated if the residual norm is less then self.getTolerance().
648    
649             :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
650             :type u0: vector value on the domain (e.g. `Data`).
651             :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
652             :type p0: scalar value on the domain (e.g. `Data`).
653             :param verbose: if set some information on iteration progress are printed
654             :type verbose: ``bool``
655             :return: flux and pressure
656             :rtype: ``tuple`` of `Data`.
657    
658             :note: The problem is solved as a least squares form
659    
660             *(I+D^*D)u+Qp=D^*f+g*
661             *Q^*u+Q^*Qp=Q^*g*
662    
663             where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
664             We eliminate the flux form the problem by setting
665    
666             *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
667    
668             form the first equation. Inserted into the second equation we get
669    
670             *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
671    
672             which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
673             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
674             """
675             self.verbose=verbose
676             rtol=self.getTolerance()
677             atol=self.getAbsoluteTolerance()
678         self.setSubProblemTolerance()
679             num_corrections=0
680             converged=False
681             p=p0
682             norm_r=None
683             while not converged:
684                   v=self.getFlux(p, fixed_flux=u0)
685                   Qp=self.__Q(p)
686                   norm_v=self.__L2(v)
687                   norm_Qp=self.__L2(Qp)
688                   if norm_v == 0.:
689                      if norm_Qp == 0.:
690                         return v,p
691                      else:
692                        fac=norm_Qp
693                   else:
694                      if norm_Qp == 0.:
695                        fac=norm_v
696                      else:
697                        fac=2./(1./norm_v+1./norm_Qp)
698                   ATOL=(atol+rtol*fac)
699                   if self.verbose:
700                        print "DarcyFlux: L2 norm of v = %e."%norm_v
701                        print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
702                        print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
703                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
704                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
705                   if norm_r == None or norm_r>ATOL:
706                       if num_corrections>max_num_corrections:
707                             raise ValueError,"maximum number of correction steps reached."
708                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
709                       num_corrections+=1
710                   else:
711                       converged=True
712             return v,p
713        def __L2(self,v):
714             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
715    
716        def __Q(self,p):
717              return util.tensor_mult(self.__permeability,util.grad(p))
718    
719        def __Aprod(self,dp):
720              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
721              Qdp=self.__Q(dp)
722              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
723              du=self.__pde_v.getSolution()
724              # self.__pde_v.getOperator().saveMM("proj.mm")
725              return Qdp+du
726        def __inner_GMRES(self,r,s):
727             return util.integrate(util.inner(r,s))
728    
729      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
730          a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
         out=-util.integrate(util.inner(a,r[0]+r[1]))  
         return out  
731    
732      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
733          if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
734          self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))            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  class StokesProblemCartesian(HomogeneousSaddlePointProblem):      def getFlux(self,p=None, fixed_flux=Data()):
739        """          """
740        Represents and solves the problem          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
741            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
742            Note that ``g`` and ``f`` are used, see `setValue`.
743    
744            :param p: pressure.
745            :type p: scalar value on the domain (e.g. `Data`).
746            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
747            :type fixed_flux: vector values on the domain (e.g. `Data`).
748            :return: flux
749            :rtype: `Data`
750            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
751                   for the permeability *k_{ij}*
752            """
753        self.setSubProblemTolerance()
754            g=self.__g
755            f=self.__f
756            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
757            if p == None:
758               self.__pde_v.setValue(Y=g)
759            else:
760               self.__pde_v.setValue(Y=g-self.__Q(p))
761            return self.__pde_v.getSolution()
762    
763        M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
764         """
765         solves
766    
767        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
768                    u_{i,i}=0
769    
770        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}            u=0 where  fixed_u_mask>0
771              eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
772    
773        If C{surface_stress} is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
774    
775        Typical usage::       typical usage:
776    
777            sp = StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
778            sp.setTolerance()              sp.setTolerance()
779            sp.initialize(...)              sp.initialize(...)
780            v,p = sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
781        """       """
782        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
783           """           """
784           Initializes the Stokes Problem.           initialize the Stokes Problem
785    
786             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
787             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
788             with macro elements for the pressure.
789    
790           @param domain: domain of the problem. The approximation order needs           :param domain: domain of the problem.
791                          to be two.           :type domain: `Domain`
          @type domain: L{Domain}  
          @warning: The approximation order needs to be two otherwise you may  
                    see oscillations in the pressure.  
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           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
803           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
804         self.__pde_proj.setValue(D=1)
805           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
806    
807        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
808           """           """
809           Assigns values to the model parameters.       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           @param f: external force       def updateStokesEquation(self, v, p):
          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar  
          @param fixed_u_mask: mask of locations with fixed velocity  
          @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}  
                              or similar  
          @param eta: viscosity  
          @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar  
          @param surface_stress: normal surface stress  
          @type surface_stress: L{Vector} object on L{FunctionSpace},  
                                L{FunctionOnBoundary} or similar  
          @param stress: initial stress  
          @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or  
                        similar  
          @note: All values need to be set.  
          """  
          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,Y=f,y=surface_stress)  
          self.__stress=stress  
   
       def B(self,v):  
          """  
          Returns M{div(v)}.  
          @return: M{div(v)}  
          @rtype: equal to the type of p  
   
          @note: Boundary conditions on p should be zero!  
          """  
          if self.show_details: print "apply divergence:"  
          self.__pde_proj.setValue(Y=-util.div(v))  
          self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
          return self.__pde_proj.getSolution(verbose=self.show_details)  
   
       def inner_pBv(self,p,Bv):  
          """  
          Returns inner product of element p and Bv (overwrite).  
   
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @return: inner product of p and Bv  
          @rtype: equal to the type of p  
855           """           """
856           s0=util.interpolate(p,Function(self.domain))           updates the Stokes equation to consider dependencies from ``v`` and ``p``
857           s1=util.interpolate(Bv,Function(self.domain))           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
858           return util.integrate(s0*s1)           """
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 inner_p(self,p0,p1):       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
891            """
892            assigns values to the model parameters
893    
894            :param f: external force
895            :type f: `Vector` object in `FunctionSpace` `Function` or similar
896            :param fixed_u_mask: mask of locations with fixed velocity.
897            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
898            :param eta: viscosity
899            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
900            :param surface_stress: normal surface stress
901            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
902            :param stress: initial stress
903        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
904            """
905            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
906    
907         def Bv(self,v,tol):
908             """
909             returns inner product of element p and div(v)
910    
911             :param v: a residual
912             :return: inner product of element p and div(v)
913             :rtype: ``float``
914             """
915             self.__pde_proj.setValue(Y=-util.div(v))
916         self.getSolverOptionsDiv().setTolerance(tol)
917         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
918             out=self.__pde_proj.getSolution()
919             return out
920    
921         def inner_pBv(self,p,Bv):
922             """
923             returns inner product of element p and Bv=-div(v)
924    
925             :param p: a pressure increment
926             :param Bv: a residual
927             :return: inner product of element p and Bv=-div(v)
928             :rtype: ``float``
929             """
930             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
931    
932         def inner_p(self,p0,p1):
933           """           """
934           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1
935    
936           @type p0: equal to the type of p           :param p0: a pressure
937           @type p1: equal to the type of p           :param p1: a pressure
938           @return: inner product of p0 and p1           :return: inner product of p0 and p1
939           @rtype: equal to the type of p           :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 inner_v(self,v0,v1):       def norm_v(self,v):
946           """           """
947           Returns inner product of two elements v0 and v1 (overwrite).           returns the norm of v
948    
949           @type v0: equal to the type of v           :param v: a velovity
950           @type v1: equal to the type of v           :return: norm of v
951           @return: inner product of v0 and v1           :rtype: non-negative ``float``
          @rtype: equal to the type of v  
952           """           """
953           gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
          gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
954    
955        def solve_A(self,u,p):  
956         def getDV(self, p, v, tol):
957           """           """
958           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).           return the value for v for a given p (overwrite)
959    
960             :param p: a pressure
961             :param v: a initial guess for the value v to return.
962             :return: dv given as *Adv=(f-Av-B^*p)*
963           """           """
964           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
965           self.__pde_u.setTolerance(self.getSubProblemTolerance())           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=-2*self.eta*util.symmetric(util.grad(u))+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-2*self.eta*util.symmetric(util.grad(u))+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):
976            """
977            Returns Bv (overwrite).
978    
979            :rtype: equal to the type of p
980            :note: boundary conditions on p should be zero!
981            """
982            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
983    
984         def solve_AinvBt(self,p, tol):
985             """
986             Solves *Av=B^*p* with accuracy `tol`
987    
988        def solve_prec(self,p):           :param p: a pressure increment
989             :return: the solution of *Av=B^*p*
990             :note: boundary conditions on v should be zero!
991           """           """
992           Applies the preconditioner.           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993             out=self.__pde_u.getSolution()
994             return  out
995    
996         def solve_prec(self,Bv, tol):
997           """           """
998           if self.show_details: print "apply preconditioner:"           applies preconditioner for for *BA^{-1}B^** to *Bv*
999           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           with accuracy `self.getSubProblemTolerance()`
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
1000    
1001             :param Bv: velocity increment
1002             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
1003             :note: boundary conditions on p are zero.
1004             """
1005             self.__pde_prec.setValue(Y=Bv)
1006         self.getSolverOptionsPressure().setTolerance(tol)
1007         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008             out=self.__pde_prec.getSolution()
1009             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26