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

revision 2349 by gross, Mon Mar 30 08:14:23 2009 UTC revision 2486 by gross, Tue Jun 23 03:38:54 2009 UTC
# 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, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38
39  class DarcyFlow(object):  class DarcyFlow(object):
# Line 48  class DarcyFlow(object): Line 48  class DarcyFlow(object):
48      @note: The problem is solved in a least squares formulation.      @note: The problem is solved in a least squares formulation.
49      """      """
50
51      def __init__(self, domain,useReduced=False):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52          """          """
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          @param domain: domain of the problem          @param domain: domain of the problem
55          @type domain: L{Domain}          @type domain: L{Domain}
56        @param useReduced: uses reduced oreder on flux and pressure
57        @type useReduced: C{bool}
58        @param adaptSubTolerance: switches on automatic subtolerance selection
60          """          """
61          self.domain=domain          self.domain=domain
62          self.__l=util.longestEdge(self.domain)**2          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          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
71          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
# Line 67  class DarcyFlow(object): Line 77  class DarcyFlow(object):
77          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78          self.setTolerance()          self.setTolerance()
79          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
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 169  class DarcyFlow(object): Line 220  class DarcyFlow(object):
220         @rtype: C{float}         @rtype: C{float}
221         """         """
222         return self.__atol         return self.__atol

def setSubProblemTolerance(self,rtol=None):
"""
Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
C{self.getTolerance()**2} is used.

@param rtol: relative tolerence
@type rtol: positive C{float}
"""
if rtol == None:
if self.getTolerance()<=0.:
raise ValueError,"A positive relative tolerance must be set."
self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
else:
if rtol<=0:
raise ValueError,"sub-problem tolerance must be positive."
self.__sub_tol=max(util.EPSILON**(0.75),rtol)

223      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
224           """      """
225           Returns the subproblem reduction factor.      Returns a suitable subtolerance
226        @type: C{float}
227           @return: subproblem reduction factor      """
228           @rtype: C{float}      return max(util.EPSILON**(0.75),self.getTolerance()**2)
229           """      def setSubProblemTolerance(self):
230           return self.__sub_tol           """
231             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
232             """
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, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
242           """           """
243           solves the problem.           solves the problem.
244
# Line 208  class DarcyFlow(object): Line 250  class DarcyFlow(object):
250           @type p0: scalar value on the domain (e.g. L{Data}).           @type p0: scalar value on the domain (e.g. L{Data}).
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
# Line 230  class DarcyFlow(object): Line 270  class DarcyFlow(object):
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 or True           self.verbose=verbose
self.show_details= show_details and self.verbose
274           rtol=self.getTolerance()           rtol=self.getTolerance()
275           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
276           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()       self.setSubProblemTolerance()
277
278           num_corrections=0           num_corrections=0
279           converged=False           converged=False
280           p=p0           p=p0
281           norm_r=None           norm_r=None
282           while not converged:           while not converged:
283                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)                 v=self.getFlux(p, fixed_flux=u0)
284                 Qp=self.__Q(p)                 Qp=self.__Q(p)
285                 norm_v=self.__L2(v)                 norm_v=self.__L2(v)
286                 norm_Qp=self.__L2(Qp)                 norm_Qp=self.__L2(Qp)
# Line 259  class DarcyFlow(object): Line 298  class DarcyFlow(object):
298                 if self.verbose:                 if self.verbose:
299                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
300                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      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                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
303                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
304                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
305                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
306                           raise ValueError,"maximum number of correction steps reached."                           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.1*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.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
308                     num_corrections+=1                     num_corrections+=1
309                 else:                 else:
310                     converged=True                     converged=True
311           return v,p           return v,p
#
#
#               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
#               #===========================================================================
#               norm_r_hat=self.__L2(r_hat)
#               norm_v=self.__L2(v)
#               norm_g=self.__L2(g)
#               norm_gv=self.__L2(g-v)
#               norm_Qp=self.__L2(Qp)
#               norm_gQp=self.__L2(g-Qp)
#               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
#               fac=min(norm_v,norm_Qp,norm_gv)
#               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
#               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
#               if r!=None:
#                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
#                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
#                   self.setSubProblemTolerance(sub_tol)
#                   print "subtol_new=",self.getSubProblemTolerance()
#               print "norm_v = ",norm_v
#               print "norm_gv = ",norm_gv
#               print "norm_Qp = ",norm_Qp
#               print "norm_gQp = ",norm_gQp
#               print "norm_g = ",norm_g
#               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
#               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
#               if fac == 0:
#                   if self.verbose: print "DarcyFlux: trivial case!"
#                   return v,p
#               #===============================================================================
#               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
#               # norm_Qp=self.__L2(Qp)
#               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
#               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
#
#               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
#               if self.verbose:
#                   print "DarcyFlux: residual = %e"%norm_r_hat
#                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
#               if norm_r_hat <= ATOL:
#                   print "DarcyFlux: iteration finalized."
#                   converged=True
#               else:
#                   # 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)
#                   # 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)
#                   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)
#               print "norm_r =",norm_r
#         return v,p
312      def __L2(self,v):      def __L2(self,v):
313           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
314
# Line 323  class DarcyFlow(object): Line 316  class DarcyFlow(object):
317
318      def __Aprod(self,dp):      def __Aprod(self,dp):
319            self.__pde_v.setTolerance(self.getSubProblemTolerance())            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
if self.show_details: print "DarcyFlux: Applying operator"
320            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
321            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
322            du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)            du=self.__pde_v.getSolution()
323              # self.__pde_v.getOperator().saveMM("proj.mm")
324            return Qdp+du            return Qdp+du
325      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
326           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 336  class DarcyFlow(object): Line 329  class DarcyFlow(object):
329           return util.integrate(util.inner(self.__Q(p), r))           return util.integrate(util.inner(self.__Q(p), r))
330
331      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
332            self.__pde_p.setTolerance(self.getSubProblemTolerance())        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
if self.show_details: print "DarcyFlux: Applying preconditioner"
333            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            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, iter_max = 100000)            # self.__pde_p.getOperator().saveMM("prec.mm")
335              return self.__pde_p.getSolution()
336
337      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):      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}          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}).          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
# Line 358  class DarcyFlow(object): Line 351  class DarcyFlow(object):
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}}          @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}}                 for the permeability M{k_{ij}}
353          """          """
354          self.__pde_v.setTolerance(self.getSubProblemTolerance())      self.setSubProblemTolerance()
355          g=self.__g          g=self.__g
356          f=self.__f          f=self.__f
357          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
# Line 366  class DarcyFlow(object): Line 359  class DarcyFlow(object):
359             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
360          else:          else:
361             self.__pde_v.setValue(Y=g-self.__Q(p))             self.__pde_v.setValue(Y=g-self.__Q(p))
362          return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)          return self.__pde_v.getSolution()
363
365       """       """
# Line 387  class StokesProblemCartesian(Homogeneous Line 380  class StokesProblemCartesian(Homogeneous
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.
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           """           """
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)
404             self.__pde_proj.setReducedOrderOn()
405         self.__pde_proj.setValue(D=1)
406             self.__pde_proj.setSymmetryOn()
407
408         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             """
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
469          """          """
470          assigns values to the model parameters          assigns values to the model parameters
# Line 423  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")
# Line 438  class StokesProblemCartesian(Homogeneous Line 494  class StokesProblemCartesian(Homogeneous
494          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
495          self.__stress=stress          self.__stress=stress
496
497       def inner_pBv(self,p,v):       def Bv(self,v):
498           """           """
499           returns inner product of element p and div(v)           returns inner product of element p and div(v)
500
# Line 447  class StokesProblemCartesian(Homogeneous Line 503  class StokesProblemCartesian(Homogeneous
503           @return: inner product of element p and div(v)           @return: inner product of element p and div(v)
504           @rtype: C{float}           @rtype: C{float}
505           """           """
506           return util.integrate(-p*util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
507             return self.__pde_proj.getSolution()
508
509         def inner_pBv(self,p,Bv):
510             """
511             returns inner product of element p and Bv=-div(v)
512
513             @param p: a pressure increment
514             @param v: a residual
515             @return: inner product of element p and Bv=-div(v)
516             @rtype: C{float}
517             """
518             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
519
520       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
521           """           """
# Line 480  class StokesProblemCartesian(Homogeneous Line 548  class StokesProblemCartesian(Homogeneous
548           @param v0: a initial guess for the value v to return.           @param v0: a initial guess for the value v to return.
549           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
550           """           """
self.__pde_u.setTolerance(self.getSubProblemTolerance())
551           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
552           if self.__stress.isEmpty():           if self.__stress.isEmpty():
553              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
554           else:           else:
555              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
556           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
557           return  out           return  out
558
559         def norm_Bv(self,Bv):
raise NotImplementedError,"no v calculation implemented."

def norm_Bv(self,v):
560          """          """
561          Returns Bv (overwrite).          Returns Bv (overwrite).
562
563          @rtype: equal to the type of p          @rtype: equal to the type of p
564          @note: boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
565          """          """
566          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
567
568       def solve_AinvBt(self,p):       def solve_AinvBt(self,p):
569           """           """
# Line 510  class StokesProblemCartesian(Homogeneous Line 573  class StokesProblemCartesian(Homogeneous
573           @return: the solution of M{Av=B^*p}           @return: the solution of M{Av=B^*p}
574           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
575           """           """
self.__pde_u.setTolerance(self.getSubProblemTolerance())
576           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
577           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
578           return  out           return  out
579
580       def solve_precB(self,v):       def solve_prec(self,Bv):
581           """           """
582           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
583           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy L{self.getSubProblemTolerance()}
584
585           @param v: velocity increment           @param v: velocity increment
586           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
587           @note: boundary conditions on p are zero.           @note: boundary conditions on p are zero.
588           """           """
589           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
590           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return self.__pde_prec.getSolution()
return self.__pde_prec.getSolution(verbose=self.show_details)

Legend:
 Removed from v.2349 changed lines Added in v.2486