/[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 2719 by gross, Wed Oct 14 06:38:03 2009 UTC revision 3515 by gross, Thu May 19 08:20:57 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 31  Some models for flow Line 32  Some models for flow
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      solves the problem     solves the problem
43      
44      *u_i+k_{ij}*p_{,j} = g_i*     *u_i+k_{ij}*p_{,j} = g_i*
45      *u_{i,i} = f*     *u_{i,i} = f*
46      
47      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49      :note: The problem is solved in a least squares formulation.     :cvar SIMPLE: simple solver
50      """     :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):     :cvar SYMSTAB: solver uses symmetric stabilization
53          """     """
54          initializes the Darcy flux problem     SIMPLE="SIMPLE"
55          :param domain: domain of the problem     POST="POST"
56          :type domain: `Domain`     STAB="STAB"
57      :param useReduced: uses reduced oreder on flux and pressure     SYMSTAB="SYMSTAB"
58      :type useReduced: ``bool``     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59      :param adaptSubTolerance: switches on automatic subtolerance selection        """
60      :type adaptSubTolerance: ``bool``          initializes the Darcy flux problem
61          """        :param domain: domain of the problem
62          self.domain=domain        :type domain: `Domain`
63          if weight == None:        :param useReduced: uses reduced oreder on flux and pressure
64             s=self.domain.getSize()        :type useReduced: ``bool``
65             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2        :param solver: solver method
66             # self.__l=(3.*util.longestEdge(self.domain))**2        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2        :param verbose: if ``True`` some information on the iteration progress are printed.
68          else:        :type verbose: ``bool``
69             self.__l=weight        :param w: weighting factor for `DarcyFlow.POST` solver
70          self.__pde_v=LinearPDESystem(domain)        :type w: ``float``
71          if useReduced: self.__pde_v.setReducedOrderOn()        
72          self.__pde_v.setSymmetryOn()        """
73          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))        self.domain=domain
74          self.__pde_p=LinearSinglePDE(domain)        self.solver=solver
75          self.__pde_p.setSymmetryOn()        self.useReduced=useReduced
76          if useReduced: self.__pde_p.setReducedOrderOn()        self.verbose=verbose
77          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        self.scale=1.
78          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        
79          self.setTolerance()        
80          self.setAbsoluteTolerance()        self.__pde_v=LinearPDESystem(domain)
81      self.__adaptSubTolerance=adaptSubTolerance        self.__pde_v.setSymmetryOn()
82      self.verbose=False        if self.useReduced: self.__pde_v.setReducedOrderOn()
83      def getSolverOptionsFlux(self):        self.__pde_p=LinearSinglePDE(domain)
84      """        self.__pde_p.setSymmetryOn()
85      Returns the solver options used to solve the flux problems        if self.useReduced: self.__pde_p.setReducedOrderOn()
86              
87      *(I+D^*D)u=F*        if self.solver  == self.SIMPLE:
88         if self.verbose: print "DarcyFlow: simple solver is used."
89             self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          elif self.solver  == self.POST:
91         self.w=w
92         if util.inf(w)<0.:
93            raise ValueError,"Weighting factor must be non-negative."
94         if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104          self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105          self.setTolerance()
106        
107            
108       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109          """
110          assigns values to model parameters
111    
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122    
123          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136          
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140            
141      :return: `SolverOptions`       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142      """           V=util.vol(self.domain)
143      return self.__pde_v.getSolverOptions()           l=V**(1./self.domain.getDim())
144      def setSolverOptionsFlux(self, options=None):          
145      """       if perm.getRank()==0:
146      Sets the solver options used to solve the flux problems          perm_inv=(1./perm)
147                    self.scale=util.integrate(perm_inv)/V*l
148      *(I+D^*D)u=F*          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149                perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150      If ``options`` is not present, the options are reset to default          
151      :param options: `SolverOptions`          
152      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.       elif perm.getRank()==2:
153      """          perm_inv=util.inverse(perm)
154      return self.__pde_v.setSolverOptions(options)              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155      def getSolverOptionsPressure(self):          perm_inv*=(1./self.scale)
156      """          perm=perm*self.scale
157      Returns the solver options used to solve the pressure problems       else:
158                raise ValueError,"illegal rank of permeability."
159      *(Q^*Q)p=Q^*G*          
160             self.__permeability=perm
161      :return: `SolverOptions`       self.__permeability_inv=perm_inv
162      """       if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163      return self.__pde_p.getSolverOptions()      
164      def setSolverOptionsPressure(self, options=None):       if self.solver  == self.SIMPLE:
165      """          self.__pde_p.setValue(A=self.__permeability)
166      Sets the solver options used to solve the pressure problems       elif self.solver  == self.POST:
167                self.__pde_p.setValue(A=self.__permeability)
168      *(Q^*Q)p=Q^*G*          k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172            self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():      
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192       def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          return self.__pde_v.getSolverOptions()
198          
199       def setSolverOptionsFlux(self, options=None):
200          """
201          Sets the solver options used to solve the flux problems
202          If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204          """
205          return self.__pde_v.setSolverOptions(options)
206        
207       def getSolverOptionsPressure(self):
208          """
209          Returns the solver options used to solve the pressure problems
210          :return: `SolverOptions`
211          """
212          return self.__pde_p.getSolverOptions()
213          
214       def setSolverOptionsPressure(self, options=None):
215          """
216          Sets the solver options used to solve the pressure problems
217          If ``options`` is not present, the options are reset to default
218          
219          :param options: `SolverOptions`
220          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221          """
222          return self.__pde_p.setSolverOptions(options)
223          
224       def setTolerance(self,rtol=1e-4):
225          """
226          sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227          
228          :param rtol: relative tolerance for the pressure
229          :type rtol: non-negative ``float``
230          """
231          if rtol<0:
232         raise ValueError,"Relative tolerance needs to be non-negative."
233          self.__rtol=rtol
234          
235       def getTolerance(self):
236          """
237          returns the relative tolerance
238          :return: current relative tolerance
239          :rtype: ``float``
240          """
241          return self.__rtol
242          
243       def solve(self,u0,p0, max_iter=100, iter_restart=20):
244          """
245          solves the problem.
246          
247          The iteration is terminated if the residual norm is less then self.getTolerance().
248    
249          :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.
250          :type u0: vector value on the domain (e.g. `escript.Data`).
251          :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.
252          :type p0: scalar value on the domain (e.g. `escript.Data`).
253          :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254          :type max_iter: ``int``
255          :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256                               A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257                               of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258          :type iter_restart: ``int``
259          :return: flux and pressure
260          :rtype: ``tuple`` of `escript.Data`.
261    
262          """
263          # rescale initial guess:
264          p0=p0/self.scale
265          if self.solver  == self.SIMPLE or self.solver  == self.POST :
266            self.__pde_p.setValue(X=self.__g ,
267                                  Y=self.__f,
268                                  y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269                                  r=p0)
270            p=self.__pde_p.getSolution()
271            u = self.getFlux(p, u0)
272          elif  self.solver  == self.STAB:
273        u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274          elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276            
277      If ``options`` is not present, the options are reset to default        if self.verbose:
278      :param options: `SolverOptions`          KGp=util.tensor_mult(self.__permeability,util.grad(p))
279      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.          def_p=self.__g-(u+KGp)
280      """          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281      return self.__pde_p.setSolverOptions(options)          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282            print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):        #rescale result
284          """        p=p*self.scale
285          assigns values to model parameters        return u,p
286          
287          :param f: volumetic sources/sinks     def getFlux(self,p, u0=None):
         :type f: scalar value on the domain (e.g. `Data`)  
         :param g: flux sources/sinks  
         :type g: vector values on the domain (e.g. `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. `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. `Data`)  
         :param permeability: permeability tensor. If scalar ``s`` is given the tensor with  
                              ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with  
                              ``v`` on the main diagonal is used.  
         :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)  
   
         :note: the values of parameters which are not set by calling ``setValue`` are not altered.  
         :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)  
                or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal  
                is along the *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 setTolerance(self,rtol=1e-4):  
         """  
         sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         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}*.  
   
         :param rtol: relative tolerance for the pressure  
         :type rtol: non-negative ``float``  
         """  
         if rtol<0:  
             raise ValueError,"Relative tolerance needs to be non-negative."  
         self.__rtol=rtol  
     def getTolerance(self):  
         """  
         returns the relative tolerance  
   
         :return: current relative tolerance  
         :rtype: ``float``  
288          """          """
289          return self.__rtol          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
   
     def setAbsoluteTolerance(self,atol=0.):  
         """  
         sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         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}*.  
   
         :param atol: absolute tolerance for the pressure  
         :type atol: non-negative ``float``  
         """  
         if atol<0:  
             raise ValueError,"Absolute tolerance needs to be non-negative."  
         self.__atol=atol  
     def getAbsoluteTolerance(self):  
        """  
        returns the absolute tolerance  
         
        :return: current absolute tolerance  
        :rtype: ``float``  
        """  
        return self.__atol  
     def getSubProblemTolerance(self):  
     """  
     Returns a suitable subtolerance  
     @type: ``float``  
     """  
     return max(util.EPSILON**(0.75),self.getTolerance()**2)  
     def setSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
          """  
      if self.__adaptSubTolerance:  
          sub_tol=self.getSubProblemTolerance()  
              self.getSolverOptionsFlux().setTolerance(sub_tol)  
          self.getSolverOptionsFlux().setAbsoluteTolerance(0.)  
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
          """  
          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 ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.  
          :type u0: vector value on the domain (e.g. `Data`).  
          :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.  
          :type p0: scalar value on the domain (e.g. `Data`).  
          :param verbose: if set some information on iteration progress are printed  
          :type verbose: ``bool``  
          :return: flux and pressure  
          :rtype: ``tuple`` of `Data`.  
   
          :note: The problem is solved as a least squares form  
   
          *(I+D^*D)u+Qp=D^*f+g*  
          *Q^*u+Q^*Qp=Q^*g*  
   
          where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
          We eliminate the flux form the problem by setting  
   
          *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  
   
          *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 *Q^*Q*). In each iteration step  
          PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
      self.setSubProblemTolerance()  
          num_corrections=0  
          converged=False  
          p=p0  
          norm_r=None  
          while not converged:  
                v=self.getFlux(p, fixed_flux=u0)  
                Qp=self.__Q(p)  
                norm_v=self.__L2(v)  
                norm_Qp=self.__L2(Qp)  
                if norm_v == 0.:  
                   if norm_Qp == 0.:  
                      return v,p  
                   else:  
                     fac=norm_Qp  
                else:  
                   if norm_Qp == 0.:  
                     fac=norm_v  
                   else:  
                     fac=2./(1./norm_v+1./norm_Qp)  
                ATOL=(atol+rtol*fac)  
                if self.verbose:  
                     print "DarcyFlux: L2 norm of v = %e."%norm_v  
                     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp  
                     print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)  
                     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)  
                     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
                if norm_r == None or norm_r>ATOL:  
                    if num_corrections>max_num_corrections:  
                          raise ValueError,"maximum number of correction steps reached."  
                    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)  
                    num_corrections+=1  
                else:  
                    converged=True  
          return v,p  
     def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
     def __Q(self,p):  
           return util.tensor_mult(self.__permeability,util.grad(p))  
   
     def __Aprod(self,dp):  
           if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
           Qdp=self.__Q(dp)  
           self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution()  
           # self.__pde_v.getOperator().saveMM("proj.mm")  
           return Qdp+du  
     def __inner_GMRES(self,r,s):  
          return util.integrate(util.inner(r,s))  
   
     def __inner_PCG(self,p,r):  
          return util.integrate(util.inner(self.__Q(p), r))  
   
     def __Msolve_PCG(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())  
           # self.__pde_p.getOperator().saveMM("prec.mm")  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=Data()):  
         """  
         returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``  
290          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` and ``f`` are used, see `setValue`.
292    
293          :param p: pressure.          :param p: pressure.
294          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
295          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          :return: flux          :return: flux
298          :rtype: `Data`          :rtype: `escript.Data`
         :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}*  
                for the permeability *k_{ij}*  
299          """          """
300      self.setSubProblemTolerance()          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301          g=self.__g              KGp=util.tensor_mult(self.__permeability,util.grad(p))
302          f=self.__f              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)              if u0 == None:
304          if p == None:             self.__pde_v.setValue(r=escript.Data())
305             self.__pde_v.setValue(Y=g)          else:
306          else:             self.__pde_v.setValue(r=u0)
307             self.__pde_v.setValue(Y=g-self.__Q(p))              u= self.__pde_v.getSolution()
308          return self.__pde_v.getSolution()      elif self.solver  == self.POST:
309                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311                if u0 == None:
312               self.__pde_v.setValue(r=escript.Data())
313            else:
314               self.__pde_v.setValue(r=u0)
315                u= self.__pde_v.getSolution()
316        elif self.solver  == self.STAB:
317             gp=util.grad(p)
318             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319                                   X= p * util.kronecker(self.domain.getDim()),
320                                   y= - p * self.domain.getNormal())                          
321             if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323             else:
324               self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326        elif  self.solver  == self.SYMSTAB:
327             gp=util.grad(p)
328             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330                                   y= escript.Data() )                          
331             if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333             else:
334               self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336        return u
337          
338        
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)  
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          p=GMRES(dp,
371                  self.__SYMSTAB_Aprod,
372              p0,
373              self.__inner,
374              atol=self.__norm(p0+dp)*self.getTolerance() ,
375              rtol=0.,
376              iter_max=max_iter,
377              iter_restart=iter_restart,
378              verbose=self.verbose,P_R=None)
379                
380              u=self.getFlux(p, u0)
381              return u,p
382    
383       def __L2(self,v):
384             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
385      
386       def __norm(self,r):
387             return util.sqrt(self.__inner(r,r))
388            
389       def __inner(self,r,s):
390             return util.integrate(util.inner(r,s), escript.Function(self.domain))
391            
392       def __STAB_Aprod(self,p):
393          gp=util.grad(p)
394          self.__pde_v.setValue(Y=-0.5*gp,
395                                X=-p*util.kronecker(self.__pde_v.getDomain()),
396                                y= p * self.domain.getNormal(),  
397                                r=escript.Data())
398          u = -self.__pde_v.getSolution()
399          self.__pde_p.setValue(Y=util.div(u),
400                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401                                y=escript.Data(),
402                                r=escript.Data())
403        
404          return  self.__pde_p.getSolution()
405      
406       def __SYMSTAB_Aprod(self,p):
407          gp=util.grad(p)
408          self.__pde_v.setValue(Y=0.5*gp ,
409                                X=escript.Data(),
410                                y=escript.Data(),  
411                                r=escript.Data())
412          u = -self.__pde_v.getSolution()
413          self.__pde_p.setValue(Y=escript.Data(),
414                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415                                y=escript.Data(),
416                                r=escript.Data())
417        
418          return  self.__pde_p.getSolution()
419          
420    
421  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
422       """       """
# Line 381  class StokesProblemCartesian(Homogeneous Line 441  class StokesProblemCartesian(Homogeneous
441           """           """
442           initialize the Stokes Problem           initialize the Stokes Problem
443    
444           :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
445             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
446             with macro elements for the pressure.
447    
448             :param domain: domain of the problem.
449           :type domain: `Domain`           :type domain: `Domain`
          :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.  
450           """           """
451           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
452           self.domain=domain           self.domain=domain
453           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
454           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
455            
456           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
457           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 405  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
468            
469       :rtype: `SolverOptions`       :rtype: `SolverOptions`
470       """       """
471       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
472       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
473           """           """
474       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 413  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476       :param options: new solver  options       :param options: new solver  options
477       :type options: `SolverOptions`       :type options: `SolverOptions`
478       """       """
479           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
480       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
481           """           """
482       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 446  class StokesProblemCartesian(Homogeneous Line 509  class StokesProblemCartesian(Homogeneous
509       """       """
510       return self.__pde_proj.getSolverOptions()       return self.__pde_proj.getSolverOptions()
511    
512       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def updateStokesEquation(self, v, p):
513             """
514             updates the Stokes equation to consider dependencies from ``v`` and ``p``
515             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
516             """
517             pass
518         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
519            """
520            assigns new values to the model parameters.
521    
522            :param f: external force
523            :type f: `Vector` object in `FunctionSpace` `Function` or similar
524            :param fixed_u_mask: mask of locations with fixed velocity.
525            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
526            :param eta: viscosity
527            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
528            :param surface_stress: normal surface stress
529            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
530            :param stress: initial stress
531        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
532            """
533            if eta !=None:
534                k=util.kronecker(self.domain.getDim())
535                kk=util.outer(k,k)
536                self.eta=util.interpolate(eta, escript.Function(self.domain))
537            self.__pde_prec.setValue(D=1/self.eta)
538                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
539            if restoration_factor!=None:
540                n=self.domain.getNormal()
541                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542            if fixed_u_mask!=None:
543                self.__pde_v.setValue(q=fixed_u_mask)
544            if f!=None: self.__f=f
545            if surface_stress!=None: self.__surface_stress=surface_stress
546            if stress!=None: self.__stress=stress
547    
548         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
549          """          """
550          assigns values to the model parameters          assigns values to the model parameters
551    
# Line 460  class StokesProblemCartesian(Homogeneous Line 559  class StokesProblemCartesian(Homogeneous
559          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
560          :param stress: initial stress          :param stress: initial stress
561      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         :note: All values needs to be set.  
562          """          """
563          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.  
         n=self.domain.getNormal()  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
564    
565       def Bv(self,v,tol):       def Bv(self,v,tol):
566           """           """
# Line 484  class StokesProblemCartesian(Homogeneous Line 570  class StokesProblemCartesian(Homogeneous
570           :return: inner product of element p and div(v)           :return: inner product of element p and div(v)
571           :rtype: ``float``           :rtype: ``float``
572           """           """
573           self.__pde_proj.setValue(Y=-util.div(v)) # -???           self.__pde_proj.setValue(Y=-util.div(v))
574       self.getSolverOptionsDiv().setTolerance(tol)       self.getSolverOptionsDiv().setTolerance(tol)
575       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
576           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
# Line 499  class StokesProblemCartesian(Homogeneous Line 585  class StokesProblemCartesian(Homogeneous
585           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
586           :rtype: ``float``           :rtype: ``float``
587           """           """
588           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
589    
590       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
591           """           """
# Line 510  class StokesProblemCartesian(Homogeneous Line 596  class StokesProblemCartesian(Homogeneous
596           :return: inner product of p0 and p1           :return: inner product of p0 and p1
597           :rtype: ``float``           :rtype: ``float``
598           """           """
599           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
600           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
601           return util.integrate(s0*s1)           return util.integrate(s0*s1)
602    
603       def norm_v(self,v):       def norm_v(self,v):
# Line 524  class StokesProblemCartesian(Homogeneous Line 610  class StokesProblemCartesian(Homogeneous
610           """           """
611           return util.sqrt(util.integrate(util.length(util.grad(v))**2))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
612    
613    
614       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
615           """           """
616           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
# Line 532  class StokesProblemCartesian(Homogeneous Line 619  class StokesProblemCartesian(Homogeneous
619           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
620           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
621           """           """
622           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.updateStokesEquation(v,p)
623             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
625       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626           if self.__stress.isEmpty():           if self.__stress.isEmpty():
627              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628           else:           else:
629              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
630           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
631           return  out           return  out
632    
633       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 549  class StokesProblemCartesian(Homogeneous Line 637  class StokesProblemCartesian(Homogeneous
637          :rtype: equal to the type of p          :rtype: equal to the type of p
638          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
639          """          """
640          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
641    
642       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
643           """           """
# Line 559  class StokesProblemCartesian(Homogeneous Line 647  class StokesProblemCartesian(Homogeneous
647           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
648           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
649           """           """
650           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
651           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
652           return  out           return  out
653    
654       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
Removed from v.2719  
changed lines
  Added in v.3515

  ViewVC Help
Powered by ViewVC 1.1.26