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

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

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

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2474 by gross, Tue Jun 16 06:32:15 2009 UTC
# Line 16  http://www.uq.edu.au/esscc Line 16  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
17  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
18  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
19  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
# Line 33  __author__="Lutz Gross, l.gross@uq.edu.a Line 33  __author__="Lutz Gross, l.gross@uq.edu.a
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      solves the problem      solves the problem
42    
43      M{u_i+k_{ij}*p_{,j} = g_i}      M{u_i+k_{ij}*p_{,j} = g_i}
44      M{u_{i,i} = f}      M{u_{i,i} = f}
45    
46      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47    
48      @note: The problem is solved in a least squares formulation.      @note: The problem is solved in a least squares formulation.
49      """      """
50    
51      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52          """          """
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          @param domain: domain of the problem          @param domain: domain of the problem
55          @type domain: L{Domain}          @type domain: L{Domain}
56        @param useReduced: uses reduced oreder on flux and pressure
57        @type useReduced: C{bool}
58        @param adaptSubTolerance: switches on automatic subtolerance selection
59        @type adaptSubTolerance: C{bool}    
60          """          """
61          self.domain=domain          self.domain=domain
62            if weight == None:
63               s=self.domain.getSize()
64               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65            else:
66               self.__l=weight
67          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
68          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
69          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
70          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
71          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
72          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
73          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
74          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
75          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
76          self.__ATOL= None          self.setTolerance()
77            self.setAbsoluteTolerance()
78        self.__adaptSubTolerance=adaptSubTolerance
79        self.verbose=False
80        def getSolverOptionsFlux(self):
81        """
82        Returns the solver options used to solve the flux problems
83        
84        M{(I+D^*D)u=F}
85        
86        @return: L{SolverOptions}
87        """
88        return self.__pde_v.getSolverOptions()
89        def setSolverOptionsFlux(self, options=None):
90        """
91        Sets the solver options used to solve the flux problems
92        
93        M{(I+D^*D)u=F}
94        
95        If C{options} is not present, the options are reset to default
96        @param options: L{SolverOptions}
97        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
98        """
99        return self.__pde_v.setSolverOptions(options)
100        def getSolverOptionsPressure(self):
101        """
102        Returns the solver options used to solve the pressure problems
103        
104        M{(Q^*Q)p=Q^*G}
105        
106        @return: L{SolverOptions}
107        """
108        return self.__pde_p.getSolverOptions()
109        def setSolverOptionsPressure(self, options=None):
110        """
111        Sets the solver options used to solve the pressure problems
112        
113        M{(Q^*Q)p=Q^*G}
114        
115        If C{options} is not present, the options are reset to default
116        @param options: L{SolverOptions}
117        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
118        """
119        return self.__pde_p.setSolverOptions(options)
120    
121      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122          """          """
# Line 78  class DarcyFlow(object): Line 130  class DarcyFlow(object):
130          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
131          @param location_of_fixed_flux:  mask for locations where flux is fixed.          @param location_of_fixed_flux:  mask for locations where flux is fixed.
132          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
133          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
134                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
135                               C{v} on the main diagonal is used.                               C{v} on the main diagonal is used.
136          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
137    
# Line 88  class DarcyFlow(object): Line 140  class DarcyFlow(object):
140                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
141                 is along the M{x_i} axis.                 is along the M{x_i} axis.
142          """          """
143          if f !=None:          if f !=None:
144             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
145             if f.isEmpty():             if f.isEmpty():
146                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
147             else:             else:
148                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
149             self.f=f             self.__f=f
150          if g !=None:            if g !=None:
151             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
152             if g.isEmpty():             if g.isEmpty():
153               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 173  class DarcyFlow(object):
173             self.__permeability=perm             self.__permeability=perm
174             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
175    
176        def setTolerance(self,rtol=1e-4):
     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}}  
177          """          """
178          self.__pde_v.setTolerance(tol)          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
179          g=self.__g  
180          f=self.__f          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
         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)  
181    
182      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
183    
184            @param rtol: relative tolerance for the pressure
185            @type rtol: non-negative C{float}
186          """          """
187          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
188          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
189          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
190                def getTolerance(self):
         @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}}  
191          """          """
192          self.__pde_v.setTolerance(tol)          returns the relative tolerance
         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)  
193    
194      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          @return: current relative tolerance
195            @rtype: C{float}
196          """          """
197          set the tolerance C{ATOL} used to terminate the solution process. It is used          return self.__rtol
   
         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  
198    
199          M{r=Q^*(g-Qp-v)}      def setAbsoluteTolerance(self,atol=0.):
200            """
201          the condition          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
202    
203          M{<(Q^*Q)^{-1} r,r> <= ATOL}          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
204    
205          holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
206    
207          @param atol: absolute tolerance for the pressure          @param atol: absolute tolerance for the pressure
208          @type atol: non-negative C{float}          @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}  
209          """          """
210          g=self.__g          if atol<0:
211          if not v_ref == None:              raise ValueError,"Absolute tolerance needs to be non-negative."
212             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)          self.__atol=atol
213          else:      def getAbsoluteTolerance(self):
214             f1=util.integrate(util.length(util.interpolate(g))**2)         """
215          if not p_ref == None:         returns the absolute tolerance
216             f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)        
217          else:         @return: current absolute tolerance
218             f2=0         @rtype: C{float}
219          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))         """
220          if self.__ATOL<=0:         return self.__atol
221             raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL      def getSubProblemTolerance(self):
222          return self.__ATOL      """
223                Returns a suitable subtolerance
224      def getTolerance(self):      @type: C{float}
225          """      """
226          returns the current tolerance.      return max(util.EPSILON**(0.75),self.getTolerance()**2)
227          def setSubProblemTolerance(self):
228          @return: used absolute tolerance.           """
229          @rtype: positive C{float}           Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
230          """           """
231          if self.__ATOL==None:       if self.__adaptSubTolerance:
232             raise ValueError,"no tolerance is defined."           sub_tol=self.getSubProblemTolerance()
233          return self.__ATOL               self.getSolverOptionsFlux().setTolerance(sub_tol)
234             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
235             self.getSolverOptionsPressure().setTolerance(sub_tol)
236             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
237             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
238    
239      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
240           """           """
241           solves the problem.           solves the problem.
242    
243           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
244    
245           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
246           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
247           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
248           @type p0: scalar value on the domain (e.g. L{Data}).           @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}  
249           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
250           @type verbose: C{bool}           @type verbose: C{bool}
          @param show_details:  if set information on the subiteration process are printed.  
          @type show_details: C{bool}  
251           @return: flux and pressure           @return: flux and pressure
252           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
253    
254           @note: The problem is solved as a least squares form           @note: The problem is solved as a least squares form
255    
256           M{(I+D^*D)u+Qp=D^*f+g}           M{(I+D^*D)u+Qp=D^*f+g}
257           M{Q^*u+Q^*Qp=Q^*g}           M{Q^*u+Q^*Qp=Q^*g}
258    
259           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
260           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
261    
262           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
263    
264           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
265    
266           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           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
267            
268           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
269           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
270           """           """
271           self.verbose=verbose           self.verbose=verbose
272           self.show_details= show_details and self.verbose           rtol=self.getTolerance()
273           self.__pde_v.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
274           self.__pde_p.setTolerance(sub_rtol)       self.setSubProblemTolerance()
275           ATOL=self.getTolerance()      
276           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL           num_corrections=0
277           #########################################################################################################################           converged=False
278           #           p=p0
279           #   we solve:           norm_r=None
280           #             while not converged:
281           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )                 v=self.getFlux(p, fixed_flux=u0)
282           #                 Qp=self.__Q(p)
283           #   residual is                 norm_v=self.__L2(v)
284           #                 norm_Qp=self.__L2(Qp)
285           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)                 if norm_v == 0.:
286           #                    if norm_Qp == 0.:
287           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                       return v,p
288           #                    else:
289           #    we use (g - Qp, v) to represent the residual. not that                      fac=norm_Qp
290           #                 else:
291           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                    if norm_Qp == 0.:
292           #                      fac=norm_v
293           #   while the initial residual is                    else:
294           #                      fac=2./(1./norm_v+1./norm_Qp)
295           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                 ATOL=(atol+rtol*fac)
296           #                   if self.verbose:
297           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))                      print "DarcyFlux: L2 norm of v = %e."%norm_v
298           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
299           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
300           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                 if norm_r == None or norm_r>ATOL:
301           self.__pde_v.setValue(r=Data())                     if num_corrections>max_num_corrections:
302           # start CG                           raise ValueError,"maximum number of correction steps reached."
303           r=ArithmeticTuple(d0, v00)                     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)
304           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)                     num_corrections+=1
305           return r[1],p                 else:
306                       converged=True
307      def __Aprod_PCG(self,dp):           return v,p
308            if self.show_details: print "DarcyFlux: Applying operator"      def __L2(self,v):
309            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
310            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))  
311            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())      def __Q(self,p):
312            du=self.__pde_v.getSolution(verbose=self.show_details)            return util.tensor_mult(self.__permeability,util.grad(p))
313            return ArithmeticTuple(mQdp,du)  
314        def __Aprod(self,dp):
315              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
316              Qdp=self.__Q(dp)
317              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
318              du=self.__pde_v.getSolution()
319              # self.__pde_v.getOperator().saveMM("proj.mm")
320              return Qdp+du
321        def __inner_GMRES(self,r,s):
322             return util.integrate(util.inner(r,s))
323    
324      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
325           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
          f0=util.integrate(util.inner(a,r[0]))  
          f1=util.integrate(util.inner(a,r[1]))  
          # print "__inner_PCG:",f0,f1,"->",f0-f1  
          return f0-f1  
326    
327      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
328            if self.show_details: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
329            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
330            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
331              return self.__pde_p.getSolution()
332    
333        def getFlux(self,p=None, fixed_flux=Data()):
334            """
335            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
336            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
337            Note that C{g} and C{f} are used, see L{setValue}.
338    
339            @param p: pressure.
340            @type p: scalar value on the domain (e.g. L{Data}).
341            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
342            @type fixed_flux: vector values on the domain (e.g. L{Data}).
343            @param tol: relative tolerance to be used.
344            @type tol: positive C{float}.
345            @return: flux
346            @rtype: L{Data}
347            @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}}
348                   for the permeability M{k_{ij}}
349            """
350        self.setSubProblemTolerance()
351            g=self.__g
352            f=self.__f
353            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
354            if p == None:
355               self.__pde_v.setValue(Y=g)
356            else:
357               self.__pde_v.setValue(Y=g-self.__Q(p))
358            return self.__pde_v.getSolution()
359    
360  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
361        """       """
362        solves       solves
363    
364            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
365                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 367  class StokesProblemCartesian(Homogeneous
367            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
368            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
369    
370        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
371    
372        typical usage:       typical usage:
373    
374              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
375              sp.setTolerance()              sp.setTolerance()
376              sp.initialize(...)              sp.initialize(...)
377              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
378        """       """
379        def __init__(self,domain,**kwargs):       def __init__(self,domain,adaptSubTolerance=True, **kwargs):
380           """           """
381           initialize the Stokes Problem           initialize the Stokes Problem
382    
383           @param domain: domain of the problem. The approximation order needs to be two.           @param domain: domain of the problem. The approximation order needs to be two.
384           @type domain: L{Domain}           @type domain: L{Domain}
385         @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
386         @type adaptSubTolerance: C{bool}
387           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
388           """           """
389           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
390           self.domain=domain           self.domain=domain
391           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
392           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
393           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
394           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
395           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
396           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
397           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
398    
399           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
400           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
401         self.__pde_proj.setValue(D=1)
402           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
403    
404        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
405             """
406         returns the solver options used  solve the equation for velocity.
407        
408         @rtype: L{SolverOptions}
409         """
410         return self.__pde_u.getSolverOptions()
411         def setSolverOptionsVelocity(self, options=None):
412             """
413         set the solver options for solving the equation for velocity.
414        
415         @param options: new solver  options
416         @type options: L{SolverOptions}
417         """
418             self.__pde_u.setSolverOptions(options)
419         def getSolverOptionsPressure(self):
420             """
421         returns the solver options used  solve the equation for pressure.
422         @rtype: L{SolverOptions}
423         """
424         return self.__pde_prec.getSolverOptions()
425         def setSolverOptionsPressure(self, options=None):
426             """
427         set the solver options for solving the equation for pressure.
428         @param options: new solver  options
429         @type options: L{SolverOptions}
430         """
431         self.__pde_prec.setSolverOptions(options)
432    
433         def setSolverOptionsDiv(self, options=None):
434             """
435         set the solver options for solving the equation to project the divergence of
436         the velocity onto the function space of presure.
437        
438         @param options: new solver options
439         @type options: L{SolverOptions}
440         """
441         self.__pde_prec.setSolverOptions(options)
442         def getSolverOptionsDiv(self):
443             """
444         returns the solver options for solving the equation to project the divergence of
445         the velocity onto the function space of presure.
446        
447         @rtype: L{SolverOptions}
448         """
449         return self.__pde_prec.getSolverOptions()
450         def setSubProblemTolerance(self):
451             """
452         Updates the tolerance for subproblems
453             """
454         if self.adaptSubTolerance():
455                 sub_tol=self.getSubProblemTolerance()
456             self.getSolverOptionsDiv().setTolerance(sub_tol)
457             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
458             self.getSolverOptionsPressure().setTolerance(sub_tol)
459             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
460             self.getSolverOptionsVelocity().setTolerance(sub_tol)
461             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
462            
463    
464         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
465          """          """
466          assigns values to the model parameters          assigns values to the model parameters
467    
# Line 383  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476          @param stress: initial stress          @param stress: initial stress
477      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
478          @note: All values needs to be set.          @note: All values needs to be set.
   
479          """          """
480          self.eta=eta          self.eta=eta
481          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
482      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
483          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
484          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
485              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
486              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
487      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
488          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
489            self.__f=f
490            self.__surface_stress=surface_stress
491          self.__stress=stress          self.__stress=stress
492    
493        def B(self,v):       def Bv(self,v):
494          """           """
495          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
496    
497          @note: boundary conditions on p should be zero!           @param p: a pressure increment
498          """           @param v: a residual
499          if self.show_details: print "apply divergence:"           @return: inner product of element p and div(v)
         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  
500           @rtype: C{float}           @rtype: C{float}
   
          @rtype: equal to the type of p  
501           """           """
502           s0=util.interpolate(p,Function(self.domain))           self.__pde_proj.setValue(Y=-util.div(v))
503           s1=util.interpolate(Bv,Function(self.domain))           return self.__pde_proj.getSolution()
          return util.integrate(s0*s1)  
504    
505        def inner_p(self,p0,p1):       def inner_pBv(self,p,Bv):
506           """           """
507           returns inner product of element p0 and p1  (overwrite)           returns inner product of element p and Bv=-div(v)
508            
509           @type p0: equal to the type of p           @param p: a pressure increment
510           @type p1: equal to the type of p           @param v: a residual
511             @return: inner product of element p and Bv=-div(v)
512           @rtype: C{float}           @rtype: C{float}
513             """
514             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
515    
516           @rtype: equal to the type of p       def inner_p(self,p0,p1):
517             """
518             Returns inner product of p0 and p1
519    
520             @param p0: a pressure
521             @param p1: a pressure
522             @return: inner product of p0 and p1
523             @rtype: C{float}
524           """           """
525           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
526           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
527           return util.integrate(s0*s1)           return util.integrate(s0*s1)
528    
529        def inner_v(self,v0,v1):       def norm_v(self,v):
530           """           """
531           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v
           
          @type v0: equal to the type of v  
          @type v1: equal to the type of v  
          @rtype: C{float}  
532    
533           @rtype: equal to the type of v           @param v: a velovity
534             @return: norm of v
535             @rtype: non-negative C{float}
536           """           """
537       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
538    
539        def solve_A(self,u,p):       def getV(self, p, v0):
540           """           """
541           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
542    
543             @param p: a pressure
544             @param v0: a initial guess for the value v to return.
545             @return: v given as M{v= A^{-1} (f-B^*p)}
546           """           """
547           if self.show_details: print "solve for velocity:"           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
548           if self.__stress.isEmpty():           if self.__stress.isEmpty():
549              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
550           else:           else:
551              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
552           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
553           return  out           return  out
554    
555        def solve_prec(self,p):       def norm_Bv(self,Bv):
556           if self.show_details: print "apply preconditioner:"          """
557           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
558           self.__pde_prec.setValue(Y=p)  
559           q=self.__pde_prec.getSolution(verbose=self.show_details)          @rtype: equal to the type of p
560           return q          @note: boundary conditions on p should be zero!
561            """
562            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
563    
564         def solve_AinvBt(self,p):
565             """
566             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
567    
568             @param p: a pressure increment
569             @return: the solution of M{Av=B^*p}
570             @note: boundary conditions on v should be zero!
571             """
572             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
573             out=self.__pde_u.getSolution()
574             return  out
575    
576         def solve_prec(self,Bv):
577             """
578             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
579             with accuracy L{self.getSubProblemTolerance()}
580    
581             @param v: velocity increment
582             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
583             @note: boundary conditions on p are zero.
584             """
585             self.__pde_prec.setValue(Y=Bv)
586             return self.__pde_prec.getSolution()

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

  ViewVC Help
Powered by ViewVC 1.1.26