/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.2208  
changed lines
  Added in v.2960

  ViewVC Help
Powered by ViewVC 1.1.26