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

revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 2351 by gross, Tue Mar 31 08:26:41 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 34  __author__="Lutz Gross, l.gross@uq.edu.a Line 34  __author__="Lutz Gross, l.gross@uq.edu.a
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
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      Represents and 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      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
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):
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          """          """
57          self.domain=domain          self.domain=domain
58            if weight == None:
59               self.__l=10.*util.longestEdge(self.domain)**2
60            else:
61               self.__l=weight
62          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
63          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          if useReduced: self.__pde_v.setReducedOrderOn()
64          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
65            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
66            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
67          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
68          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
69            if useReduced: self.__pde_p.setReducedOrderOn()
70          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
71          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
72            self.setTolerance()
73            self.setAbsoluteTolerance()
74            self.setSubProblemTolerance()
75
76      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):
77          """          """
78          Assigns values to model parameters.          assigns values to model parameters
79
80          @param f: volumetric sources/sinks          @param f: volumetic sources/sinks
81          @type f: scalar value on the domain, e.g. L{Data}          @type f: scalar value on the domain (e.g. L{Data})
82          @param g: flux sources/sinks          @param g: flux sources/sinks
83          @type g: vector value on the domain, e.g. L{Data}          @type g: vector values on the domain (e.g. L{Data})
84          @param location_of_fixed_pressure: mask for locations where pressure is fixed          @param location_of_fixed_pressure: mask for locations where pressure is fixed
85          @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})
86          @param location_of_fixed_flux: mask for locations where flux is fixed          @param location_of_fixed_flux:  mask for locations where flux is fixed.
87          @type location_of_fixed_flux: vector value on the domain (e.g. L{Data})          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
88          @param permeability: permeability tensor. If scalar C{s} is given the          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
89                               tensor with C{s} on the main diagonal is used. If                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
90                               vector C{v} is given the tensor with C{v} on the                               C{v} on the main diagonal is used.
91                               main diagonal is used.          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
92          @type permeability: scalar, vector or tensor values on the domain, e.g.
93                              L{Data}          @note: the values of parameters which are not set by calling C{setValue} are not altered.
94            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
95          @note: the values of parameters which are not set by calling                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
96                 C{setValue} are not altered                 is along the M{x_i} axis.
@note: at any point on the boundary of the domain the pressure
(C{location_of_fixed_pressure}) >0 or the normal component of
the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
the normal is along the M{x_i} axis.
97          """          """
98          if f !=None:          if f !=None:
99             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
# Line 98  class DarcyFlow(object): Line 101  class DarcyFlow(object):
101                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
102             else:             else:
103                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
104             self.f=f             self.__f=f
105          if g !=None:          if g !=None:
106             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
107             if g.isEmpty():             if g.isEmpty():
# Line 125  class DarcyFlow(object): Line 128  class DarcyFlow(object):
128             self.__permeability=perm             self.__permeability=perm
129             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))
130
131        def setTolerance(self,rtol=1e-4):
def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
132          """          """
133          Returns the flux for a given pressure C{p}.          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
134
135            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
136
137            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}}.
138
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.

@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 by
C{location_of_fixed_flux}
@type fixed_flux: vector values on the domain, e.g. L{Data}
@param tol: relative tolerance to be used
@type tol: positive float
@return: flux
@rtype: L{Data}
@note: the method uses the least squares solution
M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator
and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
"""
self.__pde_v.setTolerance(tol)
self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
return self.__pde_v.getSolution(verbose=show_details)

def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
"""
Solves the problem.

The iteration is terminated if the error in the pressure is less than
M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
side (see escript user's guide for details).

@param u0: initial guess for the flux. At locations in the domain
marked by C{location_of_fixed_flux} the value of C{u0} is
kept unchanged.
@type u0: vector value on the domain, e.g. L{Data}
@param p0: initial guess for the pressure. At locations in the domain
marked by C{location_of_fixed_pressure} the value of C{p0}
is kept unchanged.
@type p0: scalar value on the domain, e.g. L{Data}
@param atol: absolute tolerance for the pressure
@type atol: non-negative C{float}
139          @param rtol: relative tolerance for the pressure          @param rtol: relative tolerance for the pressure
140          @type rtol: non-negative C{float}          @type rtol: non-negative C{float}
141          @param sub_rtol: tolerance to be used in the sub iteration. It is          """
142                           recommended that M{sub_rtol<rtol*5.e-3}          if rtol<0:
143          @type sub_rtol: positive-negative C{float}              raise ValueError,"Relative tolerance needs to be non-negative."
144          @param verbose: if True information on iteration progress is printed          self.__rtol=rtol
145          @type verbose: C{bool}      def getTolerance(self):
146          @param show_details: if True information on the sub-iteration process          """
147                               is printed          returns the relative tolerance
148          @type show_details: C{bool}
149          @return: flux and pressure          @return: current relative tolerance
150          @rtype: C{tuple} of L{Data}          @rtype: C{float}
151            """
152          @note: the problem is solved in a least squares formulation:          return self.__rtol
153
154          M{(I+D^*D)u+Qp=D^*f+g}      def setAbsoluteTolerance(self,atol=0.):
155            """
156          M{Q^*u+Q^*Qp=Q^*g}          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if

where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the
permeability M{k_{ij}}. We eliminate the flux from the problem by
setting

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

from the first equation. Inserted into the second equation we get

M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}
on C{location_of_fixed_pressure}

which is solved using the PCG method (precondition is M{Q^*Q}).
In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}
need to be solved using a sub-iteration scheme.
"""
self.verbose=verbose
self.show_details= show_details and self.verbose
self.__pde_v.setTolerance(sub_rtol)
self.__pde_p.setTolerance(sub_rtol)
u2=u0*self.__pde_v.getCoefficient("q")
#
# first the reference velocity is calculated from
#
#   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
#
self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
u_ref=self.__pde_v.getSolution(verbose=show_details)
if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
self.__pde_v.setValue(r=Data())
#
#   and then we calculate a reference pressure
#
#       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
#
self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
p_ref=self.__pde_p.getSolution(verbose=self.show_details)
if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
self.__pde_p.setValue(r=Data())
#
#   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
#   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
#
#      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
#
#  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
#  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
#
#   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
#            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
#
#  we use the (du,Qdp) to represent the resudual
#  Q^*Q is a preconditioner
#
#  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
#
Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
if not ATOL>0:
raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)
if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
#
#   caclulate the initial residual
#
self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
du=self.__pde_v.getSolution(verbose=show_details)
r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
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)
util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
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)
157
158            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
159
160            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}}.
161
162            @param atol: absolute tolerance for the pressure
163            @type atol: non-negative C{float}
164            """
165            if atol<0:
166                raise ValueError,"Absolute tolerance needs to be non-negative."
167            self.__atol=atol
168        def getAbsoluteTolerance(self):
169           """
170           returns the absolute tolerance
171
172           @return: current absolute tolerance
173           @rtype: C{float}
174           """
175           return self.__atol
176
177        def setSubProblemTolerance(self,rtol=None):
178             """
179             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
180             C{self.getTolerance()**2} is used.
181
182             @param rtol: relative tolerence
183             @type rtol: positive C{float}
184             """
185             if rtol == None:
186                  if self.getTolerance()<=0.:
187                      raise ValueError,"A positive relative tolerance must be set."
188                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
189             else:
190                 if rtol<=0:
191                     raise ValueError,"sub-problem tolerance must be positive."
192                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
193
194        def getSubProblemTolerance(self):
195             """
196             Returns the subproblem reduction factor.
197
198             @return: subproblem reduction factor
199             @rtype: C{float}
200             """
201             return self.__sub_tol
202
203        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
204             """
205             solves the problem.
206
207             The iteration is terminated if the residual norm is less then self.getTolerance().
208
209             @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.
210             @type u0: vector value on the domain (e.g. L{Data}).
211             @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.
212             @type p0: scalar value on the domain (e.g. L{Data}).
213             @param verbose: if set some information on iteration progress are printed
214             @type verbose: C{bool}
215             @param show_details:  if set information on the subiteration process are printed.
216             @type show_details: C{bool}
217             @return: flux and pressure
218             @rtype: C{tuple} of L{Data}.
219
220             @note: The problem is solved as a least squares form
221
222             M{(I+D^*D)u+Qp=D^*f+g}
223             M{Q^*u+Q^*Qp=Q^*g}
224
225             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
226             We eliminate the flux form the problem by setting
227
228             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
229
230             form the first equation. Inserted into the second equation we get
231
232             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
233
234             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
235             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
236             """
237             self.verbose=verbose or True
238             self.show_details= show_details and self.verbose
239             rtol=self.getTolerance()
240             atol=self.getAbsoluteTolerance()
241             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
242
243             num_corrections=0
244             converged=False
245             p=p0
246             norm_r=None
247             while not converged:
248                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
249                   Qp=self.__Q(p)
250                   norm_v=self.__L2(v)
251                   norm_Qp=self.__L2(Qp)
252                   if norm_v == 0.:
253                      if norm_Qp == 0.:
254                         return v,p
255                      else:
256                        fac=norm_Qp
257                   else:
258                      if norm_Qp == 0.:
259                        fac=norm_v
260                      else:
261                        fac=2./(1./norm_v+1./norm_Qp)
262                   ATOL=(atol+rtol*fac)
263                   if self.verbose:
264                        print "DarcyFlux: L2 norm of v = %e."%norm_v
265                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
266                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
267                   if norm_r == None or norm_r>ATOL:
268                       if num_corrections>max_num_corrections:
269                             raise ValueError,"maximum number of correction steps reached."
270                       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.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
271                       num_corrections+=1
272                   else:
273                       converged=True
274             return v,p
275    #
276    #
277    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
278    #               #===========================================================================
279    #               norm_r_hat=self.__L2(r_hat)
280    #               norm_v=self.__L2(v)
281    #               norm_g=self.__L2(g)
282    #               norm_gv=self.__L2(g-v)
283    #               norm_Qp=self.__L2(Qp)
284    #               norm_gQp=self.__L2(g-Qp)
285    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
286    #               fac=min(norm_v,norm_Qp,norm_gv)
287    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
288    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
289    #               if r!=None:
290    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
291    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
292    #                   self.setSubProblemTolerance(sub_tol)
293    #                   print "subtol_new=",self.getSubProblemTolerance()
294    #               print "norm_v = ",norm_v
295    #               print "norm_gv = ",norm_gv
296    #               print "norm_Qp = ",norm_Qp
297    #               print "norm_gQp = ",norm_gQp
298    #               print "norm_g = ",norm_g
299    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
300    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
301    #               if fac == 0:
302    #                   if self.verbose: print "DarcyFlux: trivial case!"
303    #                   return v,p
304    #               #===============================================================================
305    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
306    #               # norm_Qp=self.__L2(Qp)
307    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
308    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
309    #
310    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
311    #               if self.verbose:
312    #                   print "DarcyFlux: residual = %e"%norm_r_hat
313    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
314    #               if norm_r_hat <= ATOL:
315    #                   print "DarcyFlux: iteration finalized."
316    #                   converged=True
317    #               else:
318    #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)
319    #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)
320    #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
321    #               print "norm_r =",norm_r
322    #         return v,p
323        def __L2(self,v):
324             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
325
326        def __Q(self,p):
327              return util.tensor_mult(self.__permeability,util.grad(p))
328
329        def __Aprod(self,dp):
330              self.__pde_v.setTolerance(self.getSubProblemTolerance())
331              if self.show_details: print "DarcyFlux: Applying operator"
332              Qdp=self.__Q(dp)
333              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
334              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
335              return Qdp+du
336        def __inner_GMRES(self,r,s):
337             return util.integrate(util.inner(r,s))
338
339      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
340          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
341
342      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
343          if self.show_details: print "DarcyFlux: Applying preconditioner"            self.__pde_p.setTolerance(self.getSubProblemTolerance())
344          self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))            if self.show_details: print "DarcyFlux: Applying preconditioner"
345          return self.__pde_p.getSolution(verbose=self.show_details)            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
346              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
347
348  class StokesProblemCartesian(HomogeneousSaddlePointProblem):      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
349        """          """
350        Represents and solves the problem          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
351            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
352            Note that C{g} and C{f} are used, see L{setValue}.
353
354            @param p: pressure.
355            @type p: scalar value on the domain (e.g. L{Data}).
356            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
357            @type fixed_flux: vector values on the domain (e.g. L{Data}).
358            @param tol: relative tolerance to be used.
359            @type tol: positive C{float}.
360            @return: flux
361            @rtype: L{Data}
362            @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}}
363                   for the permeability M{k_{ij}}
364            """
365            self.__pde_v.setTolerance(self.getSubProblemTolerance())
366            g=self.__g
367            f=self.__f
368            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
369            if p == None:
370               self.__pde_v.setValue(Y=g)
371            else:
372               self.__pde_v.setValue(Y=g-self.__Q(p))
373            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
374
375        M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
376         """
377         solves
378
379        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
380                    u_{i,i}=0
381
382        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}            u=0 where  fixed_u_mask>0
383              eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
384
385        If C{surface_stress} is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
386
387        Typical usage::       typical usage:
388
389            sp = StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
390            sp.setTolerance()              sp.setTolerance()
391            sp.initialize(...)              sp.initialize(...)
392            v,p = sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
393        """       """
394        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
395           """           """
396           Initializes the Stokes Problem.           initialize the Stokes Problem
397
398           @param domain: domain of the problem. The approximation order needs           @param domain: domain of the problem. The approximation order needs to be two.
to be two.
399           @type domain: L{Domain}           @type domain: L{Domain}
400           @warning: The approximation order needs to be two otherwise you may           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
see oscillations in the pressure.
401           """           """
402           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
403           self.domain=domain           self.domain=domain
# Line 311  class StokesProblemCartesian(Homogeneous Line 405  class StokesProblemCartesian(Homogeneous
405           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())
406           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
407           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
408           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
409
410           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
411           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
412           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
413           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
414
415           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
416           self.__pde_proj.setReducedOrderOn()          """
417           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
418           self.__pde_proj.setValue(D=1.)
419            @param f: external force
420        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
421           """          @param fixed_u_mask: mask of locations with fixed velocity.
422           Assigns values to the model parameters.          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
423            @param eta: viscosity
424           @param f: external force          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
425           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          @param surface_stress: normal surface stress
426           @param fixed_u_mask: mask of locations with fixed velocity          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
427           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}          @param stress: initial stress
428                               or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
429           @param eta: viscosity          @note: All values needs to be set.
430           @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar
431           @param surface_stress: normal surface stress          """
432           @type surface_stress: L{Vector} object on L{FunctionSpace},          self.eta=eta
433                                 L{FunctionOnBoundary} or similar          A =self.__pde_u.createCoefficient("A")
434           @param stress: initial stress      self.__pde_u.setValue(A=Data())
435           @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or          for i in range(self.domain.getDim()):
436                         similar          for j in range(self.domain.getDim()):
437           @note: All values need to be set.              A[i,j,j,i] += 1.
438           """              A[i,j,i,j] += 1.
439           self.eta=eta      self.__pde_prec.setValue(D=1/self.eta)
440           A =self.__pde_u.createCoefficient("A")          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
441           self.__pde_u.setValue(A=Data())          self.__f=f
442           for i in range(self.domain.getDim()):          self.__surface_stress=surface_stress
443               for j in range(self.domain.getDim()):          self.__stress=stress
444                   A[i,j,j,i] += 1.
445                   A[i,j,i,j] += 1.       def inner_pBv(self,p,v):
446           self.__pde_prec.setValue(D=1/self.eta)           """
447           self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)           returns inner product of element p and div(v)
448           self.__stress=stress
449             @param p: a pressure increment
450        def B(self,v):           @param v: a residual
451           """           @return: inner product of element p and div(v)
452           Returns M{div(v)}.           @rtype: C{float}
@return: M{div(v)}
@rtype: equal to the type of p

@note: Boundary conditions on p should be zero!
"""
if self.show_details: print "apply divergence:"
self.__pde_proj.setValue(Y=-util.div(v))
self.__pde_proj.setTolerance(self.getSubProblemTolerance())
return self.__pde_proj.getSolution(verbose=self.show_details)

def inner_pBv(self,p,Bv):
"""
Returns inner product of element p and Bv (overwrite).

@type p: equal to the type of p
@type Bv: equal to the type of result of operator B
@return: inner product of p and Bv
@rtype: equal to the type of p
453           """           """
454           s0=util.interpolate(p,Function(self.domain))           return util.integrate(-p*util.div(v))
s1=util.interpolate(Bv,Function(self.domain))
return util.integrate(s0*s1)
455
456        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
457           """           """
458           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1
459
460           @type p0: equal to the type of p           @param p0: a pressure
461           @type p1: equal to the type of p           @param p1: a pressure
462           @return: inner product of p0 and p1           @return: inner product of p0 and p1
463           @rtype: equal to the type of p           @rtype: C{float}
464           """           """
465           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
466           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
467           return util.integrate(s0*s1)           return util.integrate(s0*s1)
468
469        def inner_v(self,v0,v1):       def norm_v(self,v):
470           """           """
471           Returns inner product of two elements v0 and v1 (overwrite).           returns the norm of v
472
473           @type v0: equal to the type of v           @param v: a velovity
474           @type v1: equal to the type of v           @return: norm of v
475           @return: inner product of v0 and v1           @rtype: non-negative C{float}
@rtype: equal to the type of v
476           """           """
477           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))
478
479        def solve_A(self,u,p):       def getV(self, p, v0):
480           """           """
481           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).           return the value for v for a given p (overwrite)
482
483             @param p: a pressure
484             @param v0: a initial guess for the value v to return.
485             @return: v given as M{v= A^{-1} (f-B^*p)}
486           """           """
if self.show_details: print "solve for velocity:"
487           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
488             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
489           if self.__stress.isEmpty():           if self.__stress.isEmpty():
490              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))
491           else:           else:
492              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))
493             out=self.__pde_u.getSolution(verbose=self.show_details)
494             return  out
495
496
497             raise NotImplementedError,"no v calculation implemented."
498
499
500         def norm_Bv(self,v):
501            """
502            Returns Bv (overwrite).
503
504            @rtype: equal to the type of p
505            @note: boundary conditions on p should be zero!
506            """
507            return util.sqrt(util.integrate(util.div(v)**2))
508
509         def solve_AinvBt(self,p):
510             """
511             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
512
513             @param p: a pressure increment
514             @return: the solution of M{Av=B^*p}
515             @note: boundary conditions on v should be zero!
516             """
517             self.__pde_u.setTolerance(self.getSubProblemTolerance())
518             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
519           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
520           return out           return  out
521
522        def solve_prec(self,p):       def solve_precB(self,v):
523           """           """
524           Applies the preconditioner.           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
525             with accuracy L{self.getSubProblemTolerance()} (overwrite).
526
527             @param v: velocity increment
528             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
529             @note: boundary conditions on p are zero.
530           """           """
531           if self.show_details: print "apply preconditioner:"           self.__pde_prec.setValue(Y=-util.div(v))
532           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
533           self.__pde_prec.setValue(Y=p)           return self.__pde_prec.getSolution(verbose=self.show_details)
q=self.__pde_prec.getSolution(verbose=self.show_details)
return q

Legend:
 Removed from v.2169 changed lines Added in v.2351

 ViewVC Help Powered by ViewVC 1.1.26