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

revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2264 by gross, Wed Feb 11 06:48:28 2009 UTC
# 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      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      """      """
# Line 64  class DarcyFlow(object): Line 64  class DarcyFlow(object):
64          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
65          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67          self.__ATOL= None          self.setTolerance()
68            self.setAbsoluteTolerance()
69            self.setSubProblemTolerance()
70
71      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):
72          """          """
# Line 78  class DarcyFlow(object): Line 80  class DarcyFlow(object):
80          @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})
81          @param location_of_fixed_flux:  mask for locations where flux is fixed.          @param location_of_fixed_flux:  mask for locations where flux is fixed.
82          @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})
83          @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
84                               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
85                               C{v} on the main diagonal is used.                               C{v} on the main diagonal is used.
86          @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})
87
# Line 88  class DarcyFlow(object): Line 90  class DarcyFlow(object):
90                 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
91                 is along the M{x_i} axis.                 is along the M{x_i} axis.
92          """          """
93          if f !=None:          if f !=None:
94             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
95             if f.isEmpty():             if f.isEmpty():
96                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
97             else:             else:
98                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
99             self.f=f             self.f=f
100          if g !=None:            if g !=None:
101             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
102             if g.isEmpty():             if g.isEmpty():
103               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 123  class DarcyFlow(object):
123             self.__permeability=perm             self.__permeability=perm
124             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))
125
126        def setTolerance(self,rtol=1e-4):
def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
127          """          """
128          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
129          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
130          Note that C{g} and C{f} are used, see L{setValue}.          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
131
132          @param p: pressure.          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}}.
@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}}
"""
self.__pde_v.setTolerance(tol)
g=self.__g
f=self.__f
self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
if p == None:
self.__pde_v.setValue(Y=g)
else:
return self.__pde_v.getSolution(verbose=show_details)
133
134      def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):          @param rtol: relative tolerance for the pressure
135            @type rtol: non-negative C{float}
136          """          """
137          returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}          if rtol<0:
138          on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
139          Note that C{g} is used, see L{setValue}.          self.__rtol=rtol
140                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}}
141          """          """
142          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)
143
144      def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):          @return: current relative tolerance
145            @rtype: C{float}
146          """          """
147          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.
148
149          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      def setAbsoluteTolerance(self,atol=0.):
150            """
151          M{r=Q^*(g-Qp-v)}          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if

the condition
152
153          M{<(Q^*Q)^{-1} r,r> <= ATOL}          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
154
155          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}}.
156
157          @param atol: absolute tolerance for the pressure          @param atol: absolute tolerance for the pressure
158          @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}
159          """          """
160          g=self.__g          if atol<0:
161          if not v_ref == None:              raise ValueError,"Absolute tolerance needs to be non-negative."
162             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)          self.__atol=atol
163          else:      def getAbsoluteTolerance(self):
164             f1=util.integrate(util.length(util.interpolate(g))**2)         """
165          if not p_ref == None:         returns the absolute tolerance
167          else:         @return: current absolute tolerance
168             f2=0         @rtype: C{float}
169          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))         """
170          if self.__ATOL<=0:         return self.__atol
171             raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
172          return self.__ATOL      def setSubProblemTolerance(self,rtol=None):
173                     """
174      def getTolerance(self):           Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
175          """           C{self.getTolerance()**2} is used.
176          returns the current tolerance.
177               @param rtol: relative tolerence
178          @return: used absolute tolerance.           @type rtol: positive C{float}
179          @rtype: positive C{float}           """
180          """           if rtol == None:
181          if self.__ATOL==None:                if self.getTolerance()<=0.:
182             raise ValueError,"no tolerance is defined."                    raise ValueError,"A positive relative tolerance must be set."
183          return self.__ATOL                self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
184             else:
185                 if rtol<=0:
186                     raise ValueError,"sub-problem tolerance must be positive."
187                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
188
189      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def getSubProblemTolerance(self):
190           """           """
191             Returns the subproblem reduction factor.
192
193             @return: subproblem reduction factor
194             @rtype: C{float}
195             """
196             return self.__sub_tol
197
198        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
199             """
200           solves the problem.           solves the problem.
201
202           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().
203
204           @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.
205           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
206           @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.
207           @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}
208           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
209           @type verbose: C{bool}           @type verbose: C{bool}
210           @param show_details:  if set information on the subiteration process are printed.           @param show_details:  if set information on the subiteration process are printed.
211           @type show_details: C{bool}           @type show_details: C{bool}
212           @return: flux and pressure           @return: flux and pressure
213           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
214
215           @note: The problem is solved as a least squares form           @note: The problem is solved as a least squares form
216
217           M{(I+D^*D)u+Qp=D^*f+g}           M{(I+D^*D)u+Qp=D^*f+g}
218           M{Q^*u+Q^*Qp=Q^*g}           M{Q^*u+Q^*Qp=Q^*g}
219
220           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}}.
221           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
222
223           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
224
225           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
226
227           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
228
229           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
230           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.
231           """           """
232           self.verbose=verbose           self.verbose=verbose
233           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
234           self.__pde_v.setTolerance(sub_rtol)           rtol=self.getTolerance()
235           self.__pde_p.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
236           ATOL=self.getTolerance()           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
237           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
238           #########################################################################################################################           num_corrections=0
239           #           converged=False
240           #   we solve:           p=p0
241           #             norm_r=None
242           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )           while not converged:
243           #                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
244           #   residual is                 Qp=self.__Q(p)
245           #                 norm_v=self.__L2(v)
246           #    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)                 norm_Qp=self.__L2(Qp)
247           #                 if norm_v == 0.:
248           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC                    if norm_Qp == 0.:
249           #                       return v,p
250           #    we use (g - Qp, v) to represent the residual. not that                    else:
251           #                      fac=norm_Qp
252           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)                 else:
253           #                    if norm_Qp == 0.:
254           #   while the initial residual is                      fac=norm_v
255           #                    else:
256           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC                      fac=2./(1./norm_v+1./norm_Qp)
257           #                   ATOL=(atol+rtol*fac)
259           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)                      print "DarcyFlux: L2 norm of v = %e."%norm_v
260           v00=self.__pde_v.getSolution(verbose=show_details)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
261           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
262           self.__pde_v.setValue(r=Data())                 if norm_r == None or norm_r>ATOL:
263           # start CG                     if num_corrections>max_num_corrections:
264           r=ArithmeticTuple(d0, v00)                           raise ValueError,"maximum number of correction steps reached."
265           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)                     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)
266           return r[1],p                     num_corrections+=1
267                   else:
268                       converged=True
269             return v,p
270    #
271    #
272    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
273    #               #===========================================================================
274    #               norm_r_hat=self.__L2(r_hat)
275    #               norm_v=self.__L2(v)
276    #               norm_g=self.__L2(g)
277    #               norm_gv=self.__L2(g-v)
278    #               norm_Qp=self.__L2(Qp)
279    #               norm_gQp=self.__L2(g-Qp)
280    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
281    #               fac=min(norm_v,norm_Qp,norm_gv)
282    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
283    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
284    #               if r!=None:
285    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
286    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
287    #                   self.setSubProblemTolerance(sub_tol)
288    #                   print "subtol_new=",self.getSubProblemTolerance()
289    #               print "norm_v = ",norm_v
290    #               print "norm_gv = ",norm_gv
291    #               print "norm_Qp = ",norm_Qp
292    #               print "norm_gQp = ",norm_gQp
293    #               print "norm_g = ",norm_g
294    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
295    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
296    #               if fac == 0:
297    #                   if self.verbose: print "DarcyFlux: trivial case!"
298    #                   return v,p
299    #               #===============================================================================
300    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
301    #               # norm_Qp=self.__L2(Qp)
302    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
303    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
304    #
305    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
306    #               if self.verbose:
307    #                   print "DarcyFlux: residual = %e"%norm_r_hat
308    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
309    #               if norm_r_hat <= ATOL:
310    #                   print "DarcyFlux: iteration finalized."
311    #                   converged=True
312    #               else:
313    #                   # 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)
314    #                   # 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)
315    #                   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)
316    #               print "norm_r =",norm_r
317    #         return v,p
318        def __L2(self,v):
319             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
320
321      def __Aprod_PCG(self,dp):      def __Q(self,p):
323
324        def __Aprod(self,dp):
325              self.__pde_v.setTolerance(self.getSubProblemTolerance())
326            if self.show_details: print "DarcyFlux: Applying operator"            if self.show_details: print "DarcyFlux: Applying operator"
327            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)            Qdp=self.__Q(dp)
self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
329            du=self.__pde_v.getSolution(verbose=self.show_details)            du=self.__pde_v.getSolution(verbose=self.show_details)
330            return ArithmeticTuple(mQdp,du)            return Qdp+du
331        def __inner_GMRES(self,r,s):
332             return util.integrate(util.inner(r,s))
333
334      def __inner_PCG(self,p,r):      def __inner_PCG(self,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
336
337      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
338              self.__pde_p.setTolerance(self.getSubProblemTolerance())
339            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
340            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())
341            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
342
343
344        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
345            """
346            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
347            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
348            Note that C{g} and C{f} are used, see L{setValue}.
349
350            @param p: pressure.
351            @type p: scalar value on the domain (e.g. L{Data}).
352            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
353            @type fixed_flux: vector values on the domain (e.g. L{Data}).
354            @param tol: relative tolerance to be used.
355            @type tol: positive C{float}.
356            @return: flux
357            @rtype: L{Data}
358            @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}}
359                   for the permeability M{k_{ij}}
360            """
361            self.__pde_v.setTolerance(self.getSubProblemTolerance())
362            g=self.__g
363            f=self.__f
364            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
365            if p == None:
366               self.__pde_v.setValue(Y=g)
367            else:
368               self.__pde_v.setValue(Y=g-self.__Q(p))
369            return self.__pde_v.getSolution(verbose=show_details)
370
372        """       """
373        solves       solves
374
375            -(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}
376                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 378  class StokesProblemCartesian(Homogeneous
379            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
380
381        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
382
383        typical usage:       typical usage:
384
385              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
386              sp.setTolerance()              sp.setTolerance()
387              sp.initialize(...)              sp.initialize(...)
388              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
389        """       """
390        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
391           """           """
392           initialize the Stokes Problem           initialize the Stokes Problem
393
# Line 357  class StokesProblemCartesian(Homogeneous Line 402  class StokesProblemCartesian(Homogeneous
402           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
403           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
404           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
405
406           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
407           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
408           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
409           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
410
self.__pde_proj.setReducedOrderOn()
self.__pde_proj.setSymmetryOn()
self.__pde_proj.setValue(D=1.)

412          """          """
413          assigns values to the model parameters          assigns values to the model parameters
414
# Line 383  class StokesProblemCartesian(Homogeneous Line 423  class StokesProblemCartesian(Homogeneous
423          @param stress: initial stress          @param stress: initial stress
424      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
425          @note: All values needs to be set.          @note: All values needs to be set.
426
427          """          """
428          self.eta=eta          self.eta=eta
429          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
430      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
431          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
432          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
433              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
434              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
435      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
437            self.__f=f
438            self.__surface_stress=surface_stress
439          self.__stress=stress          self.__stress=stress
440
441        def B(self,v):       def inner_pBv(self,p,v):
442          """           """
443          returns div(v)           returns inner product of element p and div(v)
@rtype: equal to the type of p
444
445          @note: boundary conditions on p should be zero!           @param p: a pressure increment
446          """           @param v: a residual
447          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
448           @rtype: C{float}           @rtype: C{float}

@rtype: equal to the type of p
449           """           """
450           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)
451
452        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
453           """           """
454           returns inner product of element p0 and p1  (overwrite)           Returns inner product of p0 and p1

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