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

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):
        """
        initializes the Darcy flux problem
        @param domain: domain of the problem
        @type domain: L{Domain}
        """
        self.domain=domain
        self.__pde_v=LinearPDESystem(domain)
        self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
        self.__pde_v.setSymmetryOn()
        self.__pde_p=LinearSinglePDE(domain)
        self.__pde_p.setSymmetryOn()
        self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
        self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))

    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, fixed_flux=Data(),tol=1.e-8):
        """
        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, 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 float.
        @return: flux 
        @rtype: L{Data}
        @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
               for the permeability M{k_{ij}}
        """
        self.__pde_v.setTolerance(tol)
        self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)
        return self.__pde_v.getSolution()

    def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
         """ 
         solves the problem.
 
         The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where 
         C{|q|} denotes the norm of the right hand side (see escript user's guide for details).

         @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
         @type u0: vector value on the domain (e.g. L{Data}).
         @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
         @type p0: scalar value on the domain (e.g. L{Data}).
         @param atol: absolute tolerance for the pressure
         @type atol: non-negative C{float}
         @param rtol: relative tolerance for the pressure
         @type rtol: non-negative C{float}
         @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
         @type sub_rtol: positive-negative C{float}
         @param verbose: if 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)
         p2=p0*self.__pde_p.getCoefficient("q")
         u2=u0*self.__pde_v.getCoefficient("q")
         g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))
         f=self.__f-util.div(u2)
         self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())
         dv=self.__pde_v.getSolution(verbose=show_details)
         self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))
         self.__pde_p.setValue(r=Data())
         dp=self.__pde_p.getSolution(verbose=self.show_details)
         norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))
         if norm_rhs<0:
             raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."
         ATOL=util.sqrt(norm_rhs)*rtol +atol
         if not ATOL>0:
             raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)
         rhs=ArithmeticTuple(g,dv)
         dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)
         return u2+r[1],p2+dp
        
    def __Aprod_PCG(self,p):
          if self.show_details: print "DarcyFlux: Applying operator"
          Qp=util.tensor_mult(self.__permeability,util.grad(p))
          self.__pde_v.setValue(Y=Qp,X=Data())
          w=self.__pde_v.getSolution(verbose=self.show_details)
          return ArithmeticTuple(Qp,w)

    def __inner_PCG(self,p,r):
         a=util.tensor_mult(self.__permeability,util.grad(p))
         return util.integrate(util.inner(a,r[0]-r[1]))

    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]))
          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			from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG
38	gross	1414
39	gross	2100	class DarcyFlow(object):
40			"""
41			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	2100	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	2100	def __init__(self, domain):
52			"""
53			initializes the Darcy flux problem
54			@param domain: domain of the problem
55			@type domain: L{Domain}
56			"""
57			self.domain=domain
58			self.__pde_v=LinearPDESystem(domain)
59			self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
60			self.__pde_v.setSymmetryOn()
61			self.__pde_p=LinearSinglePDE(domain)
62			self.__pde_p.setSymmetryOn()
63			self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
64			self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
65	gross	1659
66	gross	2100	def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
67			"""
68			assigns values to model parameters
69	gross	1659
70	gross	2100	@param f: volumetic sources/sinks
71			@type f: scalar value on the domain (e.g. L{Data})
72			@param g: flux sources/sinks
73			@type g: vector values on the domain (e.g. L{Data})
74			@param location_of_fixed_pressure: mask for locations where pressure is fixed
75			@type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
76			@param location_of_fixed_flux: mask for locations where flux is fixed.
77			@type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
78			@param permeability: permeability tensor. If scalar C{s} is given the tensor with
79			C{s} on the main diagonal is used. If vector C{v} is given the tensor with
80			C{v} on the main diagonal is used.
81			@type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
82	gross	1659
83	gross	2100	@note: the values of parameters which are not set by calling C{setValue} are not altered.
84			@note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
85			or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
86			is along the M{x_i} axis.
87			"""
88			if f !=None:
89			f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
90			if f.isEmpty():
91			f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
92			else:
93			if f.getRank()>0: raise ValueError,"illegal rank of f."
94			self.f=f
95			if g !=None:
96			g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
97			if g.isEmpty():
98			g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
99			else:
100			if not g.getShape()==(self.domain.getDim(),):
101			raise ValueError,"illegal shape of g"
102			self.__g=g
103	gross	1659
104	gross	2100	if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
105			if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
106	gross	1659
107	gross	2100	if permeability!=None:
108			perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
109			if perm.getRank()==0:
110			perm=perm*util.kronecker(self.domain.getDim())
111			elif perm.getRank()==1:
112			perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
113			for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
114			elif perm.getRank()==2:
115			pass
116			else:
117			raise ValueError,"illegal rank of permeability."
118			self.__permeability=perm
119			self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
120	gross	1659
121
122	gross	2100	def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):
123			"""
124			returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
125			on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
126			Note that C{g} and C{f} are used, L{setValue}.
127
128			@param p: pressure.
129			@type p: scalar value on the domain (e.g. L{Data}).
130			@param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
131			@type fixed_flux: vector values on the domain (e.g. L{Data}).
132			@param tol: relative tolerance to be used.
133			@type tol: positive float.
134			@return: flux
135			@rtype: L{Data}
136			@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}}
137			for the permeability M{k_{ij}}
138			"""
139			self.__pde_v.setTolerance(tol)
140			self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)
141			return self.__pde_v.getSolution()
142	gross	1659
143	gross	2100	def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
144			"""
145			solves the problem.
146
147			The iteration is terminated if the error in the pressure is less then C{rtol * \|q\| + atol} where
148			C{\|q\|} denotes the norm of the right hand side (see escript user's guide for details).
149	gross	1659
150	gross	2100	@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.
151			@type u0: vector value on the domain (e.g. L{Data}).
152			@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.
153			@type p0: scalar value on the domain (e.g. L{Data}).
154			@param atol: absolute tolerance for the pressure
155			@type atol: non-negative C{float}
156			@param rtol: relative tolerance for the pressure
157			@type rtol: non-negative C{float}
158			@param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
159			@type sub_rtol: positive-negative C{float}
160			@param verbose: if set some information on iteration progress are printed
161			@type verbose: C{bool}
162			@param show_details: if set information on the subiteration process are printed.
163			@type show_details: C{bool}
164			@return: flux and pressure
165			@rtype: C{tuple} of L{Data}.
166
167			@note: The problem is solved as a least squares form
168	gross	1659
169	gross	2100	M{(I+D^D)u+Qp=D^f+g}
170			M{Q^u+Q^Qp=Q^*g}
171	gross	1659
172	gross	2100	where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
173			We eliminate the flux form the problem by setting
174	gross	1659
175	gross	2100	M{u=(I+D^D)^{-1}(D^f-g-Qp)} with u=u0 on location_of_fixed_flux
176
177			form the first equation. Inserted into the second equation we get
178
179			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
180
181			which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
182			PDEs with operator M{I+D^D} and with M{Q^Q} needs to be solved using a sub iteration scheme.
183			"""
184			self.verbose=verbose
185			self.show_details= show_details and self.verbose
186			self.__pde_v.setTolerance(sub_rtol)
187			self.__pde_p.setTolerance(sub_rtol)
188			p2=p0*self.__pde_p.getCoefficient("q")
189			u2=u0*self.__pde_v.getCoefficient("q")
190			g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))
191			f=self.__f-util.div(u2)
192			self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())
193			dv=self.__pde_v.getSolution(verbose=show_details)
194			self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))
195			self.__pde_p.setValue(r=Data())
196			dp=self.__pde_p.getSolution(verbose=self.show_details)
197			norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))
198			if norm_rhs<0:
199			raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."
200			ATOL=util.sqrt(norm_rhs)*rtol +atol
201			if not ATOL>0:
202			raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)
203			rhs=ArithmeticTuple(g,dv)
204	gross	2108	dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)
205	gross	2100	return u2+r[1],p2+dp
206
207			def __Aprod_PCG(self,p):
208			if self.show_details: print "DarcyFlux: Applying operator"
209			Qp=util.tensor_mult(self.__permeability,util.grad(p))
210			self.__pde_v.setValue(Y=Qp,X=Data())
211			w=self.__pde_v.getSolution(verbose=self.show_details)
212			return ArithmeticTuple(Qp,w)
213
214			def __inner_PCG(self,p,r):
215			a=util.tensor_mult(self.__permeability,util.grad(p))
216			return util.integrate(util.inner(a,r[0]-r[1]))
217
218			def __Msolve_PCG(self,r):
219			if self.show_details: print "DarcyFlux: Applying preconditioner"
220			self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]))
221			return self.__pde_p.getSolution(verbose=self.show_details)
222
223	gross	1414	class StokesProblemCartesian(HomogeneousSaddlePointProblem):
224			"""
225			solves
226
227	gross	2100	-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
228	gross	1414	u_{i,i}=0
229
230			u=0 where fixed_u_mask>0
231	gross	2100	eta(u_{i,j}+u_{j,i})n_j-p*n_i=surface_stress +stress_{ij}n_j
232	gross	1414
233	gross	2100	if surface_stress is not given 0 is assumed.
234	gross	1414
235			typical usage:
236
237			sp=StokesProblemCartesian(domain)
238			sp.setTolerance()
239			sp.initialize(...)
240			v,p=sp.solve(v0,p0)
241			"""
242			def __init__(self,domain,**kwargs):
243	gross	2100	"""
244			initialize the Stokes Problem
245
246			@param domain: domain of the problem. The approximation order needs to be two.
247			@type domain: L{Domain}
248			@warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
249			"""
250	gross	1414	HomogeneousSaddlePointProblem.__init__(self,**kwargs)
251			self.domain=domain
252			self.vol=util.integrate(1.,Function(self.domain))
253			self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
254			self.__pde_u.setSymmetryOn()
255	gross	2100	# self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
256			# self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
257	gross	1414
258			self.__pde_prec=LinearPDE(domain)
259			self.__pde_prec.setReducedOrderOn()
260	gross	2100	self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
261	gross	1414	self.__pde_prec.setSymmetryOn()
262
263			self.__pde_proj=LinearPDE(domain)
264			self.__pde_proj.setReducedOrderOn()
265			self.__pde_proj.setSymmetryOn()
266			self.__pde_proj.setValue(D=1.)
267
268	gross	2100	def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
269			"""
270			assigns values to the model parameters
271
272			@param f: external force
273			@type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
274			@param fixed_u_mask: mask of locations with fixed velocity.
275			@type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
276			@param eta: viscosity
277			@type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
278			@param surface_stress: normal surface stress
279			@type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
280			@param stress: initial stress
281			@type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
282			@note: All values needs to be set.
283
284			"""
285	gross	1414	self.eta=eta
286	gross	1841	A =self.__pde_u.createCoefficient("A")
287	gross	1414	self.__pde_u.setValue(A=Data())
288			for i in range(self.domain.getDim()):
289			for j in range(self.domain.getDim()):
290			A[i,j,j,i] += 1.
291			A[i,j,i,j] += 1.
292	artak	1554	self.__pde_prec.setValue(D=1/self.eta)
293	gross	1414	self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
294	gross	2100	self.__stress=stress
295	gross	1414
296	gross	2100	def B(self,v):
297			"""
298			returns div(v)
299			@rtype: equal to the type of p
300	gross	1414
301	gross	2100	@note: boundary conditions on p should be zero!
302			"""
303			if self.show_details: print "apply divergence:"
304			self.__pde_proj.setValue(Y=-util.div(v))
305			self.__pde_proj.setTolerance(self.getSubProblemTolerance())
306			return self.__pde_proj.getSolution(verbose=self.show_details)
307
308			def inner_pBv(self,p,Bv):
309			"""
310			returns inner product of element p and Bv (overwrite)
311
312			@type p: equal to the type of p
313			@type Bv: equal to the type of result of operator B
314			@rtype: C{float}
315
316			@rtype: equal to the type of p
317			"""
318			s0=util.interpolate(p,Function(self.domain))
319			s1=util.interpolate(Bv,Function(self.domain))
320	gross	1414	return util.integrate(s0*s1)
321
322	gross	2100	def inner_p(self,p0,p1):
323			"""
324			returns inner product of element p0 and p1 (overwrite)
325
326			@type p0: equal to the type of p
327			@type p1: equal to the type of p
328			@rtype: C{float}
329	artak	1550
330	gross	2100	@rtype: equal to the type of p
331			"""
332			s0=util.interpolate(p0/self.eta,Function(self.domain))
333			s1=util.interpolate(p1/self.eta,Function(self.domain))
334			return util.integrate(s0*s1)
335	artak	1550
336	gross	2100	def inner_v(self,v0,v1):
337			"""
338			returns inner product of two element v0 and v1 (overwrite)
339
340			@type v0: equal to the type of v
341			@type v1: equal to the type of v
342			@rtype: C{float}
343	gross	1414
344	gross	2100	@rtype: equal to the type of v
345			"""
346			gv0=util.grad(v0)
347			gv1=util.grad(v1)
348			return util.integrate(util.inner(gv0,gv1))
349
350	gross	1414	def solve_A(self,u,p):
351			"""
352			solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
353			"""
354	gross	2100	if self.show_details: print "solve for velocity:"
355	gross	1414	self.__pde_u.setTolerance(self.getSubProblemTolerance())
356	gross	2100	if self.__stress.isEmpty():
357			self.__pde_u.setValue(X=-2self.etautil.symmetric(util.grad(u))+p*util.kronecker(self.domain))
358			else:
359			self.__pde_u.setValue(X=self.__stress-2self.etautil.symmetric(util.grad(u))+p*util.kronecker(self.domain))
360			out=self.__pde_u.getSolution(verbose=self.show_details)
361			return out
362	gross	1414
363			def solve_prec(self,p):
364	gross	2100	if self.show_details: print "apply preconditioner:"
365	gross	1414	self.__pde_prec.setTolerance(self.getSubProblemTolerance())
366			self.__pde_prec.setValue(Y=p)
367			q=self.__pde_prec.getSolution(verbose=self.show_details)
368	artak	1554	return q