/[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 2445 by gross, Fri May 29 03:23:25 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, weight=None, 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          if weight == None:          if weight == None:
# Line 64  class DarcyFlow(object): Line 68  class DarcyFlow(object):
68          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
69          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
70          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
         # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)  
71          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
72          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
73          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
# Line 72  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 174  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        """
226        return max(util.EPSILON**(0.75),self.getTolerance()**2)
227        def setSubProblemTolerance(self):
228             """
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           @return: subproblem reduction factor      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
          @rtype: C{float}  
          """  
          return self.__sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):  
240           """           """
241           solves the problem.           solves the problem.
242    
# Line 213  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 236  class DarcyFlow(object): Line 269  class DarcyFlow(object):
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           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 280  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, iter_max = 100000)            du=self.__pde_v.getSolution()
319            # self.__pde_v.getOperator().saveMM("proj.mm")            # 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):
# Line 294  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            # self.__pde_p.getOperator().saveMM("prec.mm")            # self.__pde_p.getOperator().saveMM("prec.mm")
331            return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)            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 317  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 325  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, iter_max=100000)          return self.__pde_v.getSolution()
359    
360  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
361       """       """
# Line 346  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.ILU0)  
   
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)           self.__pde_proj=LinearPDE(domain)
# Line 372  class StokesProblemCartesian(Homogeneous Line 401  class StokesProblemCartesian(Homogeneous
401       self.__pde_proj.setValue(D=1)       self.__pde_proj.setValue(D=1)
402           self.__pde_proj.setSymmetryOn()           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          """          """
# Line 388  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 413  class StokesProblemCartesian(Homogeneous Line 500  class StokesProblemCartesian(Homogeneous
500           @rtype: C{float}           @rtype: C{float}
501           """           """
502           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
503           return self.__pde_proj.getSolution()           return self.__pde_proj.getSolution()
504    
505       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
# Line 458  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):       def norm_Bv(self,Bv):
# Line 484  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_prec(self,Bv):       def solve_prec(self,Bv):
# Line 499  class StokesProblemCartesian(Homogeneous Line 583  class StokesProblemCartesian(Homogeneous
583           @note: boundary conditions on p are zero.           @note: boundary conditions on p are zero.
584           """           """
585           self.__pde_prec.setValue(Y=Bv)           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)  
   
 # rename solve_prec and change argument v to Bv  
 # chnage the argument of inner_pBv to v->Bv  
 # add Bv  
 # inner p still needed?  
 # change norm_Bv argument to Bv  

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

  ViewVC Help
Powered by ViewVC 1.1.26