/[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 2349 by gross, Mon Mar 30 08:14:23 2009 UTC revision 2620 by gross, Thu Aug 20 06:24:00 2009 UTC
# Line 1  Line 1 
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
5  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
6  #  #
# Line 10  Line 10 
10  #  #
11  ########################################################  ########################################################
12    
13  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
15  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# 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               # 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()
80          self.setSubProblemTolerance()      self.__adaptSubTolerance=adaptSubTolerance
81        self.verbose=False
82        def getSolverOptionsFlux(self):
83        """
84        Returns the solver options used to solve the flux problems
85        
86        M{(I+D^*D)u=F}
87        
88        @return: L{SolverOptions}
89        """
90        return self.__pde_v.getSolverOptions()
91        def setSolverOptionsFlux(self, options=None):
92        """
93        Sets the solver options used to solve the flux problems
94        
95        M{(I+D^*D)u=F}
96        
97        If C{options} is not present, the options are reset to default
98        @param options: L{SolverOptions}
99        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
100        """
101        return self.__pde_v.setSolverOptions(options)
102        def getSolverOptionsPressure(self):
103        """
104        Returns the solver options used to solve the pressure problems
105        
106        M{(Q^*Q)p=Q^*G}
107        
108        @return: L{SolverOptions}
109        """
110        return self.__pde_p.getSolverOptions()
111        def setSolverOptionsPressure(self, options=None):
112        """
113        Sets the solver options used to solve the pressure problems
114        
115        M{(Q^*Q)p=Q^*G}
116        
117        If C{options} is not present, the options are reset to default
118        @param options: L{SolverOptions}
119        @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
120        """
121        return self.__pde_p.setSolverOptions(options)
122    
123      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124          """          """
# Line 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             """
233         if self.__adaptSubTolerance:
234             sub_tol=self.getSubProblemTolerance()
235                 self.getSolverOptionsFlux().setTolerance(sub_tol)
236             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
237             self.getSolverOptionsPressure().setTolerance(sub_tol)
238             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
239             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
240    
241      def solve(self,u0,p0, 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):
316            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
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    
364  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
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.
390         @type adaptSubTolerance: C{bool}
391           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
392           """           """
393           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
394           self.domain=domain           self.domain=domain
395           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
396           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
397           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
398           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
399           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
400           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
401           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
402    
403       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):           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             """
458         if self.adaptSubTolerance():
459                 sub_tol=self.getSubProblemTolerance()
460             self.getSolverOptionsDiv().setTolerance(sub_tol)
461             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
462             self.getSolverOptionsPressure().setTolerance(sub_tol)
463             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
464             self.getSolverOptionsVelocity().setTolerance(sub_tol)
465             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
466            
467    
468         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
469          """          """
470          assigns values to the model parameters          assigns values to the model parameters
471    
# Line 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 432  class StokesProblemCartesian(Homogeneous Line 488  class StokesProblemCartesian(Homogeneous
488          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
489              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
490              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
491            n=self.domain.getNormal()
492      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
493          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
494          self.__f=f          self.__f=f
495          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
496          self.__stress=stress          self.__stress=stress
497    
498       def inner_pBv(self,p,v):       def Bv(self,v):
499           """           """
500           returns inner product of element p and div(v)           returns inner product of element p and div(v)
501    
# Line 447  class StokesProblemCartesian(Homogeneous Line 504  class StokesProblemCartesian(Homogeneous
504           @return: inner product of element p and div(v)           @return: inner product of element p and div(v)
505           @rtype: C{float}           @rtype: C{float}
506           """           """
507           return util.integrate(-p*util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
508             return self.__pde_proj.getSolution()
509    
510         def inner_pBv(self,p,Bv):
511             """
512             returns inner product of element p and Bv=-div(v)
513    
514             @param p: a pressure increment
515             @param v: a residual
516             @return: inner product of element p and Bv=-div(v)
517             @rtype: C{float}
518             """
519             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
520    
521       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
522           """           """
# Line 480  class StokesProblemCartesian(Homogeneous Line 549  class StokesProblemCartesian(Homogeneous
549           @param v0: a initial guess for the value v to return.           @param v0: a initial guess for the value v to return.
550           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
551           """           """
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
552           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)
553           if self.__stress.isEmpty():           if self.__stress.isEmpty():
554              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
555           else:           else:
556              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
557           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
558           return  out           return  out
559    
560         def norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
561          """          """
562          Returns Bv (overwrite).          Returns Bv (overwrite).
563    
564          @rtype: equal to the type of p          @rtype: equal to the type of p
565          @note: boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
566          """          """
567          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
568    
569       def solve_AinvBt(self,p):       def solve_AinvBt(self,p):
570           """           """
# Line 510  class StokesProblemCartesian(Homogeneous Line 574  class StokesProblemCartesian(Homogeneous
574           @return: the solution of M{Av=B^*p}           @return: the solution of M{Av=B^*p}
575           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
576           """           """
          self.__pde_u.setTolerance(self.getSubProblemTolerance())  
577           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))
578           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution()
579           return  out           return  out
580    
581       def solve_precB(self,v):       def solve_prec(self,Bv):
582           """           """
583           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
584           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy L{self.getSubProblemTolerance()}
585    
586           @param v: velocity increment           @param v: velocity increment
587           @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^*}
588           @note: boundary conditions on p are zero.           @note: boundary conditions on p are zero.
589           """           """
590           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
591           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.2620

  ViewVC Help
Powered by ViewVC 1.1.26