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

Legend:
Removed from v.2156  
changed lines
  Added in v.2620

  ViewVC Help
Powered by ViewVC 1.1.26