escript/py_src/flows.py

########################################################
#
# Copyright (c) 2003-2008 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2003-2008 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="http://www.uq.edu.au/esscc/escript-finley"

"""
Some models for flow

@var __author__: name of author
@var __copyright__: copyrights
@var __license__: licence agreement
@var __url__: url entry point on documentation
@var __version__: version
@var __date__: date of the version
"""

__author__="Lutz Gross, l.gross@uq.edu.au"

from escript import *
import util
from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm

class DarcyFlow(object):
    """
    solves the problem 

    M{u_i+k_{ij}*p_{,j} = g_i}
    M{u_{i,i} = f}

    where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability, 

    @note: The problem is solved in a least squares formulation.
    """

    def __init__(self, domain,useReduced=False):
        """
        initializes the Darcy flux problem
        @param domain: domain of the problem
        @type domain: L{Domain}
        """
        self.domain=domain
        self.__pde_v=LinearPDESystem(domain)
        if useReduced: self.__pde_v.setReducedOrderOn()
        self.__pde_v.setSymmetryOn()
        self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
        self.__pde_p=LinearSinglePDE(domain)
        self.__pde_p.setSymmetryOn()
        if useReduced: self.__pde_p.setReducedOrderOn()
        self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
        self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
        self.__ATOL= None

    def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
        """
        assigns values to model parameters

        @param f: volumetic sources/sinks
        @type f: scalar value on the domain (e.g. L{Data})
        @param g: flux sources/sinks
        @type g: vector values on the domain (e.g. L{Data})
        @param location_of_fixed_pressure: mask for locations where pressure is fixed
        @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
        @param location_of_fixed_flux:  mask for locations where flux is fixed.
        @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
        @param permeability: permeability tensor. If scalar C{s} is given the tensor with 
                             C{s} on the main diagonal is used. If vector C{v} is given the tensor with 
                             C{v} on the main diagonal is used.
        @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})

        @note: the values of parameters which are not set by calling C{setValue} are not altered.
        @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 direction of the normal
               is along the M{x_i} axis.
        """
        if f !=None: 
           f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
           if f.isEmpty():
               f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
           else:
               if f.getRank()>0: raise ValueError,"illegal rank of f."
           self.f=f
        if g !=None:   
           g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
           if g.isEmpty():
             g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
           else:
             if not g.getShape()==(self.domain.getDim(),):
               raise ValueError,"illegal shape of g"
           self.__g=g

        if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
        if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)

        if permeability!=None:
           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
           if perm.getRank()==0:
               perm=perm*util.kronecker(self.domain.getDim())
           elif perm.getRank()==1:
               perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
               for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
           elif perm.getRank()==2:
              pass
           else:
              raise ValueError,"illegal rank of permeability."
           self.__permeability=perm
           self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))


    def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
        """
        returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
        on locations where C{location_of_fixed_flux} is positive (see L{setValue}). 
        Note that C{g} and C{f} are used, see L{setValue}.
        
        @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}}
        """
        self.__pde_v.setTolerance(tol)
        g=self.__g
        f=self.__f
        self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
        if p == None:
           self.__pde_v.setValue(Y=g)
        else:
           self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
        return self.__pde_v.getSolution(verbose=show_details)

    def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
        """
        returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
        on locations where C{location_of_fixed_pressure} is positive (see L{setValue}). 
        Note that C{g} is used, see L{setValue}.
        
        @param v: flux.
        @type v: vector-valued on the domain (e.g. L{Data}).
        @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
        @type fixed_pressure: vector values on the domain (e.g. L{Data}).
        @param tol: relative tolerance to be used.
        @type tol: positive C{float}.
        @return: pressure 
        @rtype: L{Data}
        @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}}
               for the permeability M{k_{ij}}
        """
        self.__pde_v.setTolerance(tol)
        g=self.__g
        self.__pde_p.setValue(r=fixed_pressure)
        if v == None:
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
        else:
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
        return self.__pde_p.getSolution(verbose=show_details)

    def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
        """
        set the tolerance C{ATOL} used to terminate the solution process. It is used

        M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}

        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.

        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

        M{r=Q^*(g-Qp-v)}

        the condition 

        M{<(Q^*Q)^{-1} r,r> <= ATOL} 

        holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}

        @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 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}
        """
        g=self.__g
        if not v_ref == None:
           f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
        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
        
    def getTolerance(self):
        """
        returns the current tolerance.
   
        @return: used absolute tolerance.
        @rtype: positive C{float}
        """
        if self.__ATOL==None:
           raise ValueError,"no tolerance is defined."
        return self.__ATOL

    def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
         """ 
         solves the problem.
 
         The iteration is terminated if the residual norm is less then self.getTolerance().

         @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 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 set some information on iteration progress are printed
         @type verbose: C{bool}
         @param show_details:  if set information on the subiteration process are printed.
         @type show_details: C{bool}
         @return: flux and pressure
         @rtype: C{tuple} of L{Data}.
 
         @note: The problem is solved as a least squares form

         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 form the problem by setting 

         M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux

         form 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 p=p0  on 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} needs 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)
         ATOL=self.getTolerance()
         if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
         #########################################################################################################################
         # 
         #   we solve:
         #  
         #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) ) 
         #
         #   residual is 
         #
         #    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)
         #
         #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC 
         #
         #    we use (g - Qp, v) to represent the residual. not that 
         #
         #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
         #
         #   while the initial residual is 
         #
         #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
         #   
         d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
         self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
         v00=self.__pde_v.getSolution(verbose=show_details)
         if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
         self.__pde_v.setValue(r=Data())
         # start CG
         r=ArithmeticTuple(d0, v00)
         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)
         return r[1],p

    def __Aprod_PCG(self,dp):
          if self.show_details: print "DarcyFlux: Applying operator"
          #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
          mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
          self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
          du=self.__pde_v.getSolution(verbose=self.show_details)
          return ArithmeticTuple(mQdp,du)

    def __inner_PCG(self,p,r):
         a=util.tensor_mult(self.__permeability,util.grad(p))
         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

    def __Msolve_PCG(self,r):
          if self.show_details: print "DarcyFlux: Applying preconditioner"
          self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
          return self.__pde_p.getSolution(verbose=self.show_details)

class StokesProblemCartesian(HomogeneousSaddlePointProblem):
      """
      solves 

          -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
                u_{i,i}=0

          u=0 where  fixed_u_mask>0
          eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j

      if surface_stress is not given 0 is assumed. 

      typical usage:

            sp=StokesProblemCartesian(domain)
            sp.setTolerance()
            sp.initialize(...)
            v,p=sp.solve(v0,p0)
      """
      def __init__(self,domain,**kwargs):
         """
         initialize the Stokes Problem

         @param domain: domain of the problem. The approximation order needs to be two.
         @type domain: L{Domain}
         @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
         """
         HomogeneousSaddlePointProblem.__init__(self,**kwargs)
         self.domain=domain
         self.vol=util.integrate(1.,Function(self.domain))
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()
         # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
            
         self.__pde_prec=LinearPDE(domain)
         self.__pde_prec.setReducedOrderOn()
         # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
         self.__pde_prec.setSymmetryOn()

         self.__pde_proj=LinearPDE(domain)
         self.__pde_proj.setReducedOrderOn()
         self.__pde_proj.setSymmetryOn()
         self.__pde_proj.setValue(D=1.)

      def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
        """
        assigns values to the model parameters

        @param f: external force
        @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 eta: 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 needs 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 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}

         @rtype: equal to the type of p
         """
         s0=util.interpolate(p,Function(self.domain))
         s1=util.interpolate(Bv,Function(self.domain))
         return util.integrate(s0*s1)

      def inner_p(self,p0,p1):
         """
         returns inner product of element p0 and p1  (overwrite)
         
         @type p0: equal to the type of p
         @type p1: equal to the type of p
         @rtype: C{float}

         @rtype: equal to the type of p
         """
         s0=util.interpolate(p0/self.eta,Function(self.domain))
         s1=util.interpolate(p1/self.eta,Function(self.domain))
         return util.integrate(s0*s1)

      def inner_v(self,v0,v1):
         """
         returns inner product of two element v0 and v1  (overwrite)
         
         @type v0: equal to the type of v
         @type v1: equal to the type of v
         @rtype: C{float}

         @rtype: equal to the type of v
         """
     gv0=util.grad(v0)
     gv1=util.grad(v1)
         return util.integrate(util.inner(gv0,gv1))

      def solve_A(self,u,p):
         """
         solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
         """
         if self.show_details: print "solve for velocity:"
         self.__pde_u.setTolerance(self.getSubProblemTolerance())
         if self.__stress.isEmpty():
            self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
         else:
            self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
         out=self.__pde_u.getSolution(verbose=self.show_details)
         return  out

      def solve_prec(self,p):
         if self.show_details: print "apply preconditioner:"
         self.__pde_prec.setTolerance(self.getSubProblemTolerance())
         self.__pde_prec.setValue(Y=p)
         q=self.__pde_prec.getSolution(verbose=self.show_details)
         return q
1	ksteube	1809	########################################################
2	gross	1414	#
3	ksteube	1809	# Copyright (c) 2003-2008 by University of Queensland
4			# Earth Systems Science Computational Center (ESSCC)
5			# http://www.uq.edu.au/esscc
6	gross	1414	#
7	ksteube	1809	# Primary Business: Queensland, Australia
8			# Licensed under the Open Software License version 3.0
9			# http://www.opensource.org/licenses/osl-3.0.php
10	gross	1414	#
11	ksteube	1809	########################################################
12	gross	1414
13	ksteube	1809	__copyright__="""Copyright (c) 2003-2008 by University of Queensland
14			Earth Systems Science Computational Center (ESSCC)
15			http://www.uq.edu.au/esscc
16			Primary Business: Queensland, Australia"""
17			__license__="""Licensed under the Open Software License version 3.0
18			http://www.opensource.org/licenses/osl-3.0.php"""
19			__url__="http://www.uq.edu.au/esscc/escript-finley"
20
21	gross	1414	"""
22			Some models for flow
23
24			@var __author__: name of author
25			@var __copyright__: copyrights
26			@var __license__: licence agreement
27			@var __url__: url entry point on documentation
28			@var __version__: version
29			@var __date__: date of the version
30			"""
31
32			__author__="Lutz Gross, l.gross@uq.edu.au"
33
34			from escript import *
35			import util
36	gross	2100	from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37	gross	2156	from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
38	gross	1414
39	gross	2100	class DarcyFlow(object):
40			"""
41	gross	2208	solves the problem
42	gross	1659
43	gross	2100	M{u_i+k_{ij}*p_{,j} = g_i}
44			M{u_{i,i} = f}
45	gross	1659
46	gross	2208	where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47	gross	1659
48	gross	2100	@note: The problem is solved in a least squares formulation.
49			"""
50	gross	1659
51	gross	2208	def __init__(self, domain,useReduced=False):
52	gross	2100	"""
53	gross	2208	initializes the Darcy flux problem
54	gross	2100	@param domain: domain of the problem
55			@type domain: L{Domain}
56			"""
57			self.domain=domain
58			self.__pde_v=LinearPDESystem(domain)
59	gross	2208	if useReduced: self.__pde_v.setReducedOrderOn()
60			self.__pde_v.setSymmetryOn()
61	gross	2100	self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
62			self.__pde_p=LinearSinglePDE(domain)
63			self.__pde_p.setSymmetryOn()
64	gross	2208	if useReduced: self.__pde_p.setReducedOrderOn()
65	gross	2100	self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66			self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67	gross	2208	self.__ATOL= None
68	gross	1659
69	gross	2100	def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70			"""
71	gross	2208	assigns values to model parameters
72	gross	1659
73	gross	2208	@param f: volumetic sources/sinks
74			@type f: scalar value on the domain (e.g. L{Data})
75	gross	2100	@param g: flux sources/sinks
76	gross	2208	@type g: vector values on the domain (e.g. L{Data})
77	gross	2100	@param location_of_fixed_pressure: mask for locations where pressure is fixed
78	gross	2208	@type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
79			@param location_of_fixed_flux: mask for locations where flux is fixed.
80			@type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
81			@param permeability: permeability tensor. If scalar C{s} is given the tensor with
82			C{s} on the main diagonal is used. If vector C{v} is given the tensor with
83			C{v} on the main diagonal is used.
84			@type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
85	gross	1659
86	gross	2208	@note: the values of parameters which are not set by calling C{setValue} are not altered.
87			@note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
88			or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
89			is along the M{x_i} axis.
90	gross	2100	"""
91	gross	2208	if f !=None:
92	gross	2100	f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
93			if f.isEmpty():
94			f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
95			else:
96			if f.getRank()>0: raise ValueError,"illegal rank of f."
97			self.f=f
98	gross	2208	if g !=None:
99	gross	2100	g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
100			if g.isEmpty():
101			g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
102			else:
103			if not g.getShape()==(self.domain.getDim(),):
104			raise ValueError,"illegal shape of g"
105			self.__g=g
106	gross	1659
107	gross	2100	if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
108			if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
109	gross	1659
110	gross	2100	if permeability!=None:
111			perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
112			if perm.getRank()==0:
113			perm=perm*util.kronecker(self.domain.getDim())
114			elif perm.getRank()==1:
115			perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
116			for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
117			elif perm.getRank()==2:
118			pass
119			else:
120			raise ValueError,"illegal rank of permeability."
121			self.__permeability=perm
122			self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
123	gross	1659
124
125	gross	2208	def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
126	gross	2100	"""
127	gross	2208	returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
128			on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
129			Note that C{g} and C{f} are used, see L{setValue}.
130
131			@param p: pressure.
132			@type p: scalar value on the domain (e.g. L{Data}).
133			@param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134			@type fixed_flux: vector values on the domain (e.g. L{Data}).
135			@param tol: relative tolerance to be used.
136			@type tol: positive C{float}.
137			@return: flux
138	gross	2100	@rtype: L{Data}
139	gross	2208	@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}}
140			for the permeability M{k_{ij}}
141	gross	2100	"""
142			self.__pde_v.setTolerance(tol)
143	gross	2208	g=self.__g
144			f=self.__f
145			self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
146			if p == None:
147			self.__pde_v.setValue(Y=g)
148			else:
149			self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
150	gross	2156	return self.__pde_v.getSolution(verbose=show_details)
151	gross	1659
152	gross	2208	def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
153	caltinay	2169	"""
154	gross	2208	returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
155			on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).
156			Note that C{g} is used, see L{setValue}.
157
158			@param v: flux.
159			@type v: vector-valued on the domain (e.g. L{Data}).
160			@param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
161			@type fixed_pressure: vector values on the domain (e.g. L{Data}).
162			@param tol: relative tolerance to be used.
163			@type tol: positive C{float}.
164			@return: pressure
165			@rtype: L{Data}
166			@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}}
167			for the permeability M{k_{ij}}
168			"""
169			self.__pde_v.setTolerance(tol)
170			g=self.__g
171			self.__pde_p.setValue(r=fixed_pressure)
172			if v == None:
173			self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
174			else:
175			self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
176			return self.__pde_p.getSolution(verbose=show_details)
177	gross	1659
178	gross	2208	def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
179			"""
180			set the tolerance C{ATOL} used to terminate the solution process. It is used
181	gross	1659
182	gross	2208	M{ATOL = atol + rtol * max( \|g-v_ref\|, \|Qp_ref\| )}
183
184			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.
185
186			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
187
188			M{r=Q^*(g-Qp-v)}
189
190			the condition
191
192			M{<(Q^*Q)^{-1} r,r> <= ATOL}
193
194			holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
195
196	caltinay	2169	@param atol: absolute tolerance for the pressure
197			@type atol: non-negative C{float}
198			@param rtol: relative tolerance for the pressure
199			@type rtol: non-negative C{float}
200	gross	2208	@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
201			L{getPressure} method or use the value from a previous time step.
202			@type p_ref: scalar value on the domain (e.g. L{Data}).
203			@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
204			L{getFlux} method or use the value from a previous time step.
205			@type v_ref: vector-valued on the domain (e.g. L{Data}).
206			@return: used absolute tolerance.
207			@rtype: positive C{float}
208			"""
209			g=self.__g
210			if not v_ref == None:
211			f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212			else:
213			f1=util.integrate(util.length(util.interpolate(g))**2)
214			if not p_ref == None:
215			f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
216			else:
217			f2=0
218			self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219			if self.__ATOL<=0:
220			raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221			return self.__ATOL
222
223			def getTolerance(self):
224			"""
225			returns the current tolerance.
226
227			@return: used absolute tolerance.
228			@rtype: positive C{float}
229			"""
230			if self.__ATOL==None:
231			raise ValueError,"no tolerance is defined."
232			return self.__ATOL
233	gross	1659
234	gross	2208	def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235			"""
236			solves the problem.
237
238			The iteration is terminated if the residual norm is less then self.getTolerance().
239	gross	1659
240	gross	2208	@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.
241			@type u0: vector value on the domain (e.g. L{Data}).
242			@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.
243			@type p0: scalar value on the domain (e.g. L{Data}).
244			@param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245			@type sub_rtol: positive-negative C{float}
246			@param verbose: if set some information on iteration progress are printed
247			@type verbose: C{bool}
248			@param show_details: if set information on the subiteration process are printed.
249			@type show_details: C{bool}
250			@return: flux and pressure
251			@rtype: C{tuple} of L{Data}.
252
253			@note: The problem is solved as a least squares form
254	gross	2100
255	gross	2208	M{(I+D^D)u+Qp=D^f+g}
256			M{Q^u+Q^Qp=Q^*g}
257	gross	2100
258	gross	2208	where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
259			We eliminate the flux form the problem by setting
260	caltinay	2169
261	gross	2208	M{u=(I+D^D)^{-1}(D^f-g-Qp)} with u=u0 on location_of_fixed_flux
262	caltinay	2169
263	gross	2208	form the first equation. Inserted into the second equation we get
264	caltinay	2169
265	gross	2208	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
266
267			which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
268			PDEs with operator M{I+D^D} and with M{Q^Q} needs to be solved using a sub iteration scheme.
269			"""
270			self.verbose=verbose
271			self.show_details= show_details and self.verbose
272			self.__pde_v.setTolerance(sub_rtol)
273			self.__pde_p.setTolerance(sub_rtol)
274			ATOL=self.getTolerance()
275			if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276			#########################################################################################################################
277			#
278			# we solve:
279			#
280			# Q^(I-(I+D^D)^{-1})Q dp = Q^* (g-u0-Qp0 - (I+D^D)^{-1} ( D^(f-Du0)+g-u0-Qp0) )
281			#
282			# residual is
283			#
284			# 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)
285			#
286			# with v = (I+D^D)^{-1} (D^f+g-Qp) including BC
287			#
288			# we use (g - Qp, v) to represent the residual. not that
289			#
290			# dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291			#
292			# while the initial residual is
293			#
294			# r0=( g - Qp0, v00) with v00=(I+D^D)^{-1} (D^f+g-Qp0) including BC
295			#
296			d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297			self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
298			v00=self.__pde_v.getSolution(verbose=show_details)
299			if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300			self.__pde_v.setValue(r=Data())
301			# start CG
302			r=ArithmeticTuple(d0, v00)
303			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)
304			return r[1],p
305	caltinay	2169
306	gross	2208	def __Aprod_PCG(self,dp):
307			if self.show_details: print "DarcyFlux: Applying operator"
308			# -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
309			mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
310			self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311			du=self.__pde_v.getSolution(verbose=self.show_details)
312			return ArithmeticTuple(mQdp,du)
313	caltinay	2169
314	gross	2100	def __inner_PCG(self,p,r):
315	gross	2208	a=util.tensor_mult(self.__permeability,util.grad(p))
316			f0=util.integrate(util.inner(a,r[0]))
317			f1=util.integrate(util.inner(a,r[1]))
318			# print "__inner_PCG:",f0,f1,"->",f0-f1
319			return f0-f1
320	gross	2100
321			def __Msolve_PCG(self,r):
322	gross	2208	if self.show_details: print "DarcyFlux: Applying preconditioner"
323			self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
324			return self.__pde_p.getSolution(verbose=self.show_details)
325	gross	2100
326	gross	1414	class StokesProblemCartesian(HomogeneousSaddlePointProblem):
327			"""
328	gross	2208	solves
329	gross	1414
330	gross	2208	-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
331			u_{i,i}=0
332	gross	1414
333	gross	2208	u=0 where fixed_u_mask>0
334			eta(u_{i,j}+u_{j,i})n_j-p*n_i=surface_stress +stress_{ij}n_j
335	gross	1414
336	gross	2208	if surface_stress is not given 0 is assumed.
337	gross	1414
338	gross	2208	typical usage:
339	gross	1414
340	gross	2208	sp=StokesProblemCartesian(domain)
341			sp.setTolerance()
342			sp.initialize(...)
343			v,p=sp.solve(v0,p0)
344	gross	1414	"""
345			def __init__(self,domain,**kwargs):
346	gross	2100	"""
347	gross	2208	initialize the Stokes Problem
348	gross	2100
349	gross	2208	@param domain: domain of the problem. The approximation order needs to be two.
350	gross	2100	@type domain: L{Domain}
351	gross	2208	@warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
352	gross	2100	"""
353	gross	1414	HomogeneousSaddlePointProblem.__init__(self,**kwargs)
354			self.domain=domain
355			self.vol=util.integrate(1.,Function(self.domain))
356			self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
357			self.__pde_u.setSymmetryOn()
358	gross	2100	# self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
359			# self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
360	gross	2208
361	gross	1414	self.__pde_prec=LinearPDE(domain)
362			self.__pde_prec.setReducedOrderOn()
363	gross	2156	# self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364	gross	1414	self.__pde_prec.setSymmetryOn()
365
366			self.__pde_proj=LinearPDE(domain)
367			self.__pde_proj.setReducedOrderOn()
368			self.__pde_proj.setSymmetryOn()
369			self.__pde_proj.setValue(D=1.)
370
371	gross	2100	def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
372	gross	2208	"""
373			assigns values to the model parameters
374	gross	2100
375	gross	2208	@param f: external force
376			@type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
377			@param fixed_u_mask: mask of locations with fixed velocity.
378			@type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
379			@param eta: viscosity
380			@type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
381			@param surface_stress: normal surface stress
382			@type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
383			@param stress: initial stress
384			@type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
385			@note: All values needs to be set.
386
387			"""
388			self.eta=eta
389			A =self.__pde_u.createCoefficient("A")
390			self.__pde_u.setValue(A=Data())
391			for i in range(self.domain.getDim()):
392			for j in range(self.domain.getDim()):
393			A[i,j,j,i] += 1.
394			A[i,j,i,j] += 1.
395			self.__pde_prec.setValue(D=1/self.eta)
396			self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
397			self.__stress=stress
398	gross	1414
399	gross	2100	def B(self,v):
400	gross	2208	"""
401			returns div(v)
402			@rtype: equal to the type of p
403	gross	1414
404	gross	2208	@note: boundary conditions on p should be zero!
405			"""
406			if self.show_details: print "apply divergence:"
407			self.__pde_proj.setValue(Y=-util.div(v))
408			self.__pde_proj.setTolerance(self.getSubProblemTolerance())
409			return self.__pde_proj.getSolution(verbose=self.show_details)
410	gross	2100
411			def inner_pBv(self,p,Bv):
412			"""
413	gross	2208	returns inner product of element p and Bv (overwrite)
414
415	gross	2100	@type p: equal to the type of p
416			@type Bv: equal to the type of result of operator B
417	gross	2208	@rtype: C{float}
418
419	gross	2100	@rtype: equal to the type of p
420			"""
421			s0=util.interpolate(p,Function(self.domain))
422			s1=util.interpolate(Bv,Function(self.domain))
423	gross	1414	return util.integrate(s0*s1)
424
425	gross	2100	def inner_p(self,p0,p1):
426			"""
427	gross	2208	returns inner product of element p0 and p1 (overwrite)
428
429	gross	2100	@type p0: equal to the type of p
430			@type p1: equal to the type of p
431	gross	2208	@rtype: C{float}
432
433	gross	2100	@rtype: equal to the type of p
434			"""
435			s0=util.interpolate(p0/self.eta,Function(self.domain))
436			s1=util.interpolate(p1/self.eta,Function(self.domain))
437			return util.integrate(s0*s1)
438	artak	1550
439	gross	2100	def inner_v(self,v0,v1):
440			"""
441	gross	2208	returns inner product of two element v0 and v1 (overwrite)
442
443	gross	2100	@type v0: equal to the type of v
444			@type v1: equal to the type of v
445	gross	2208	@rtype: C{float}
446
447	gross	2100	@rtype: equal to the type of v
448			"""
449	gross	2208	gv0=util.grad(v0)
450			gv1=util.grad(v1)
451	gross	2100	return util.integrate(util.inner(gv0,gv1))
452
453	gross	1414	def solve_A(self,u,p):
454			"""
455	gross	2208	solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
456	gross	1414	"""
457	gross	2100	if self.show_details: print "solve for velocity:"
458	gross	1414	self.__pde_u.setTolerance(self.getSubProblemTolerance())
459	gross	2100	if self.__stress.isEmpty():
460			self.__pde_u.setValue(X=-2self.etautil.symmetric(util.grad(u))+p*util.kronecker(self.domain))
461			else:
462			self.__pde_u.setValue(X=self.__stress-2self.etautil.symmetric(util.grad(u))+p*util.kronecker(self.domain))
463			out=self.__pde_u.getSolution(verbose=self.show_details)
464	gross	2208	return out
465	gross	1414
466			def solve_prec(self,p):
467	gross	2100	if self.show_details: print "apply preconditioner:"
468	gross	1414	self.__pde_prec.setTolerance(self.getSubProblemTolerance())
469			self.__pde_prec.setValue(Y=p)
470			q=self.__pde_prec.getSolution(verbose=self.show_details)
471	artak	1554	return q