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

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

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

revision 2344 by jfenwick, Mon Mar 30 02:13:58 2009 UTC revision 2474 by gross, Tue Jun 16 06:32:15 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
59        @type adaptSubTolerance: C{bool}    
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            else:
66               self.__l=weight
67          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
68          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
69          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
# Line 67  class DarcyFlow(object): Line 75  class DarcyFlow(object):
75          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
76          self.setTolerance()          self.setTolerance()
77          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
78          self.setSubProblemTolerance()      self.__adaptSubTolerance=adaptSubTolerance
79        self.verbose=False
80        def getSolverOptionsFlux(self):
81        """
82        Returns the solver options used to solve the flux problems
83        
84        M{(I+D^*D)u=F}
85        
86        @return: L{SolverOptions}
87        """
88        return self.__pde_v.getSolverOptions()
89        def setSolverOptionsFlux(self, options=None):
90        """
91        Sets the solver options used to solve the flux problems
92        
93        M{(I+D^*D)u=F}
94        
95        If C{options} is not present, the options are reset to default
96        @param options: L{SolverOptions}
97        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
98        """
99        return self.__pde_v.setSolverOptions(options)
100        def getSolverOptionsPressure(self):
101        """
102        Returns the solver options used to solve the pressure problems
103        
104        M{(Q^*Q)p=Q^*G}
105        
106        @return: L{SolverOptions}
107        """
108        return self.__pde_p.getSolverOptions()
109        def setSolverOptionsPressure(self, options=None):
110        """
111        Sets the solver options used to solve the pressure problems
112        
113        M{(Q^*Q)p=Q^*G}
114        
115        If C{options} is not present, the options are reset to default
116        @param options: L{SolverOptions}
117        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
118        """
119        return self.__pde_p.setSolverOptions(options)
120    
121      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
122          """          """
# Line 169  class DarcyFlow(object): Line 218  class DarcyFlow(object):
218         @rtype: C{float}         @rtype: C{float}
219         """         """
220         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)  
   
221      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
222           """      """
223           Returns the subproblem reduction factor.      Returns a suitable subtolerance
224        @type: C{float}
225           @return: subproblem reduction factor      """
226           @rtype: C{float}      return max(util.EPSILON**(0.75),self.getTolerance()**2)
227           """      def setSubProblemTolerance(self):
228           return self.__sub_tol           """
229             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
230             """
231         if self.__adaptSubTolerance:
232             sub_tol=self.getSubProblemTolerance()
233                 self.getSolverOptionsFlux().setTolerance(sub_tol)
234             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
235             self.getSolverOptionsPressure().setTolerance(sub_tol)
236             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
237             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
238    
239      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
240           """           """
241           solves the problem.           solves the problem.
242    
# Line 208  class DarcyFlow(object): Line 248  class DarcyFlow(object):
248           @type p0: scalar value on the domain (e.g. L{Data}).           @type p0: scalar value on the domain (e.g. L{Data}).
249           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
250           @type verbose: C{bool}           @type verbose: C{bool}
          @param show_details:  if set information on the subiteration process are printed.  
          @type show_details: C{bool}  
251           @return: flux and pressure           @return: flux and pressure
252           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
253    
# Line 230  class DarcyFlow(object): Line 268  class DarcyFlow(object):
268           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
269           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
270           """           """
271           self.verbose=verbose or True           self.verbose=verbose
          self.show_details= show_details and self.verbose  
272           rtol=self.getTolerance()           rtol=self.getTolerance()
273           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
274           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()       self.setSubProblemTolerance()
275        
276           num_corrections=0           num_corrections=0
277           converged=False           converged=False
278           p=p0           p=p0
279           norm_r=None           norm_r=None
280           while not converged:           while not converged:
281                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)                 v=self.getFlux(p, fixed_flux=u0)
282                 Qp=self.__Q(p)                 Qp=self.__Q(p)
283                 norm_v=self.__L2(v)                 norm_v=self.__L2(v)
284                 norm_Qp=self.__L2(Qp)                 norm_Qp=self.__L2(Qp)
# Line 263  class DarcyFlow(object): Line 300  class DarcyFlow(object):
300                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
301                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
302                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
303                     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)
304                     num_corrections+=1                     num_corrections+=1
305                 else:                 else:
306                     converged=True                     converged=True
307           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  
308      def __L2(self,v):      def __L2(self,v):
309           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))
310    
# Line 323  class DarcyFlow(object): Line 312  class DarcyFlow(object):
312            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
313    
314      def __Aprod(self,dp):      def __Aprod(self,dp):
315            self.__pde_v.setTolerance(self.getSubProblemTolerance())            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
           if self.show_details: print "DarcyFlux: Applying operator"  
316            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
317            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
318            du=self.__pde_v.getSolution(verbose=self.show_details)            du=self.__pde_v.getSolution()
319              # self.__pde_v.getOperator().saveMM("proj.mm")
320            return Qdp+du            return Qdp+du
321      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
322           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 336  class DarcyFlow(object): Line 325  class DarcyFlow(object):
325           return util.integrate(util.inner(self.__Q(p), r))           return util.integrate(util.inner(self.__Q(p), r))
326    
327      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
328            self.__pde_p.setTolerance(self.getSubProblemTolerance())        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
           if self.show_details: print "DarcyFlux: Applying preconditioner"  
329            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())
330            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
331              return self.__pde_p.getSolution()
332    
333      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):      def getFlux(self,p=None, fixed_flux=Data()):
334          """          """
335          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
336          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
# Line 358  class DarcyFlow(object): Line 347  class DarcyFlow(object):
347          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
348                 for the permeability M{k_{ij}}                 for the permeability M{k_{ij}}
349          """          """
350          self.__pde_v.setTolerance(self.getSubProblemTolerance())      self.setSubProblemTolerance()
351          g=self.__g          g=self.__g
352          f=self.__f          f=self.__f
353          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 355  class DarcyFlow(object):
355             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
356          else:          else:
357             self.__pde_v.setValue(Y=g-self.__Q(p))             self.__pde_v.setValue(Y=g-self.__Q(p))
358          return self.__pde_v.getSolution(verbose=show_details)          return self.__pde_v.getSolution()
359    
360  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
361       """       """
# Line 387  class StokesProblemCartesian(Homogeneous Line 376  class StokesProblemCartesian(Homogeneous
376              sp.initialize(...)              sp.initialize(...)
377              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
378       """       """
379       def __init__(self,domain,**kwargs):       def __init__(self,domain,adaptSubTolerance=True, **kwargs):
380           """           """
381           initialize the Stokes Problem           initialize the Stokes Problem
382    
383           @param domain: domain of the problem. The approximation order needs to be two.           @param domain: domain of the problem. The approximation order needs to be two.
384           @type domain: L{Domain}           @type domain: L{Domain}
385         @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
386         @type adaptSubTolerance: C{bool}
387           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
388           """           """
389           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
390           self.domain=domain           self.domain=domain
391           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
392           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
393           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
394           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
395           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
396           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
397           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
398    
399             self.__pde_proj=LinearPDE(domain)
400             self.__pde_proj.setReducedOrderOn()
401         self.__pde_proj.setValue(D=1)
402             self.__pde_proj.setSymmetryOn()
403    
404         def getSolverOptionsVelocity(self):
405             """
406         returns the solver options used  solve the equation for velocity.
407        
408         @rtype: L{SolverOptions}
409         """
410         return self.__pde_u.getSolverOptions()
411         def setSolverOptionsVelocity(self, options=None):
412             """
413         set the solver options for solving the equation for velocity.
414        
415         @param options: new solver  options
416         @type options: L{SolverOptions}
417         """
418             self.__pde_u.setSolverOptions(options)
419         def getSolverOptionsPressure(self):
420             """
421         returns the solver options used  solve the equation for pressure.
422         @rtype: L{SolverOptions}
423         """
424         return self.__pde_prec.getSolverOptions()
425         def setSolverOptionsPressure(self, options=None):
426             """
427         set the solver options for solving the equation for pressure.
428         @param options: new solver  options
429         @type options: L{SolverOptions}
430         """
431         self.__pde_prec.setSolverOptions(options)
432    
433         def setSolverOptionsDiv(self, options=None):
434             """
435         set the solver options for solving the equation to project the divergence of
436         the velocity onto the function space of presure.
437        
438         @param options: new solver options
439         @type options: L{SolverOptions}
440         """
441         self.__pde_prec.setSolverOptions(options)
442         def getSolverOptionsDiv(self):
443             """
444         returns the solver options for solving the equation to project the divergence of
445         the velocity onto the function space of presure.
446        
447         @rtype: L{SolverOptions}
448         """
449         return self.__pde_prec.getSolverOptions()
450         def setSubProblemTolerance(self):
451             """
452         Updates the tolerance for subproblems
453             """
454         if self.adaptSubTolerance():
455                 sub_tol=self.getSubProblemTolerance()
456             self.getSolverOptionsDiv().setTolerance(sub_tol)
457             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
458             self.getSolverOptionsPressure().setTolerance(sub_tol)
459             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
460             self.getSolverOptionsVelocity().setTolerance(sub_tol)
461             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
462            
463    
464       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
465          """          """
466          assigns values to the model parameters          assigns values to the model parameters
# Line 423  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476          @param stress: initial stress          @param stress: initial stress
477      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
478          @note: All values needs to be set.          @note: All values needs to be set.
   
479          """          """
480          self.eta=eta          self.eta=eta
481          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
# Line 438  class StokesProblemCartesian(Homogeneous Line 490  class StokesProblemCartesian(Homogeneous
490          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
491          self.__stress=stress          self.__stress=stress
492    
493       def inner_pBv(self,p,v):       def Bv(self,v):
494           """           """
495           returns inner product of element p and div(v)           returns inner product of element p and div(v)
496    
# Line 447  class StokesProblemCartesian(Homogeneous Line 499  class StokesProblemCartesian(Homogeneous
499           @return: inner product of element p and div(v)           @return: inner product of element p and div(v)
500           @rtype: C{float}           @rtype: C{float}
501           """           """
502           return util.integrate(-p*util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
503             return self.__pde_proj.getSolution()
504    
505         def inner_pBv(self,p,Bv):
506             """
507             returns inner product of element p and Bv=-div(v)
508    
509             @param p: a pressure increment
510             @param v: a residual
511             @return: inner product of element p and Bv=-div(v)
512             @rtype: C{float}
513             """
514             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
515    
516       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
517           """           """
# Line 480  class StokesProblemCartesian(Homogeneous Line 544  class StokesProblemCartesian(Homogeneous
544           @param v0: a initial guess for the value v to return.           @param v0: a initial guess for the value v to return.
545           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
546           """           """
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
547           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)
548           if self.__stress.isEmpty():           if self.__stress.isEmpty():
549              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
550           else:           else:
551              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
552           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
553           return  out           return  out
554    
555         def norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
556          """          """
557          Returns Bv (overwrite).          Returns Bv (overwrite).
558    
559          @rtype: equal to the type of p          @rtype: equal to the type of p
560          @note: boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
561          """          """
562          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
563    
564       def solve_AinvBt(self,p):       def solve_AinvBt(self,p):
565           """           """
# Line 510  class StokesProblemCartesian(Homogeneous Line 569  class StokesProblemCartesian(Homogeneous
569           @return: the solution of M{Av=B^*p}           @return: the solution of M{Av=B^*p}
570           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
571           """           """
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
572           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))
573           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
574           return  out           return  out
575    
576       def solve_precB(self,v):       def solve_prec(self,Bv):
577           """           """
578           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
579           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy L{self.getSubProblemTolerance()}
580    
581           @param v: velocity increment           @param v: velocity increment
582           @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^*}
583           @note: boundary conditions on p are zero.           @note: boundary conditions on p are zero.
584           """           """
585           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
586           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return self.__pde_prec.getSolution()
          return self.__pde_prec.getSolution(verbose=self.show_details)  

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

  ViewVC Help
Powered by ViewVC 1.1.26