/[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 3360 by jfenwick, Thu Nov 18 00:20:21 2010 UTC revision 3432 by jfenwick, Fri Jan 7 01:32:07 2011 UTC
# Line 32  Some models for flow Line 32  Some models for flow
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
# Line 80  class DarcyFlow(object): Line 80  class DarcyFlow(object):
80        # if self.useReduced: self.__pde_l.setReducedOrderOn()        # if self.useReduced: self.__pde_l.setReducedOrderOn()
81        self.setTolerance()        self.setTolerance()
82        self.setAbsoluteTolerance()        self.setAbsoluteTolerance()
83        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85                
86     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
87        """        """
# Line 176  class DarcyFlow(object): Line 176  class DarcyFlow(object):
176        if f !=None:        if f !=None:
177       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178       if f.isEmpty():       if f.isEmpty():
179          f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))          f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180       else:       else:
181          if f.getRank()>0: raise ValueError,"illegal rank of f."          if f.getRank()>0: raise ValueError,"illegal rank of f."
182          self.__f=f          self.__f=f
183        if g !=None:        if g !=None:
184       g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))       g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185       if g.isEmpty():       if g.isEmpty():
186          g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))          g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187       else:       else:
188          if not g.getShape()==(self.domain.getDim(),):          if not g.getShape()==(self.domain.getDim(),):
189             raise ValueError,"illegal shape of g"             raise ValueError,"illegal shape of g"
# Line 241  class DarcyFlow(object): Line 241  class DarcyFlow(object):
241    
242     def __Aprod_v(self,dv):     def __Aprod_v(self,dv):
243        # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
244        dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)        dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245        a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp        a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
246        b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)        b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247        return ArithmeticTuple(a,b,-dp)        return ArithmeticTuple(a,b,-dp)
# Line 249  class DarcyFlow(object): Line 249  class DarcyFlow(object):
249     def __Msolve_PCG_v(self,r):     def __Msolve_PCG_v(self,r):
250        # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv        # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251        if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252        self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())        self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
253        # self.__pde_p.getOperator().saveMM("prec.mm")        # self.__pde_p.getOperator().saveMM("prec.mm")
254        return self.__pde_k.getSolution()        return self.__pde_k.getSolution()
255    
# Line 259  class DarcyFlow(object): Line 259  class DarcyFlow(object):
259     def __Aprod_p(self,dp):     def __Aprod_p(self,dp):
260        if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"        if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261        Gdp=util.grad(dp)        Gdp=util.grad(dp)
262        self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())        self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
263        du=self.__pde_k.getSolution()        du=self.__pde_k.getSolution()
264        # self.__pde_v.getOperator().saveMM("proj.mm")        # self.__pde_v.getOperator().saveMM("proj.mm")
265        return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)        return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
# Line 279  class DarcyFlow(object): Line 279  class DarcyFlow(object):
279        #v=self.__getFlux(p, u0, f=self.__f, g=g2)              #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
280     def __Msolve_PCG_p(self,r):     def __Msolve_PCG_p(self,r):
281        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282        self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())        self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
283        # self.__pde_p.getOperator().saveMM("prec.mm")        # self.__pde_p.getOperator().saveMM("prec.mm")
284        return self.__pde_p.getSolution()        return self.__pde_p.getSolution()
285                    
# Line 287  class DarcyFlow(object): Line 287  class DarcyFlow(object):
287         return util.integrate(util.inner(util.grad(p), r[0]-r[1]))         return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288    
289     def __L2(self,v):     def __L2(self,v):
290        return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))        return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
291    
292     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293        """        """
# Line 474  class DarcyFlowOld(object): Line 474  class DarcyFlowOld(object):
474          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
475          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
476          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
477          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
478          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
479          self.setTolerance()          self.setTolerance()
480          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
481      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 546  class DarcyFlowOld(object): Line 546  class DarcyFlowOld(object):
546          if f !=None:          if f !=None:
547             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
548             if f.isEmpty():             if f.isEmpty():
549                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
550             else:             else:
551                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
552             self.__f=f             self.__f=f
553          if g !=None:          if g !=None:
554             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
555             if g.isEmpty():             if g.isEmpty():
556               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
557             else:             else:
558               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
559                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 698  class DarcyFlowOld(object): Line 698  class DarcyFlowOld(object):
698                 if self.verbose:                 if self.verbose:
699                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
700                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
701                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
702                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
703                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
704                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
705                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
706                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
707                     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)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.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)
708                     num_corrections+=1                     num_corrections+=1
709                 else:                 else:
710                     converged=True                     converged=True
711           return v,p           return v,p
712      def __L2(self,v):      def __L2(self,v):
713           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
714    
715      def __Q(self,p):      def __Q(self,p):
716            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 718  class DarcyFlowOld(object): Line 718  class DarcyFlowOld(object):
718      def __Aprod(self,dp):      def __Aprod(self,dp):
719            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
720            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
721            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
722            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
723            # self.__pde_v.getOperator().saveMM("proj.mm")            # self.__pde_v.getOperator().saveMM("proj.mm")
724            return Qdp+du            return Qdp+du
# Line 730  class DarcyFlowOld(object): Line 730  class DarcyFlowOld(object):
730    
731      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
732        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
733            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=escript.Data(), r=escript.Data())
734            # self.__pde_p.getOperator().saveMM("prec.mm")            # self.__pde_p.getOperator().saveMM("prec.mm")
735            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
736    
737      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
738          """          """
739          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
740          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
# Line 874  class StokesProblemCartesian(Homogeneous Line 874  class StokesProblemCartesian(Homogeneous
874          if eta !=None:          if eta !=None:
875              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
876              kk=util.outer(k,k)              kk=util.outer(k,k)
877              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
878          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
879              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
880          if restoration_factor!=None:          if restoration_factor!=None:
# Line 886  class StokesProblemCartesian(Homogeneous Line 886  class StokesProblemCartesian(Homogeneous
886          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
887          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
888    
889       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
890          """          """
891          assigns values to the model parameters          assigns values to the model parameters
892    
# Line 926  class StokesProblemCartesian(Homogeneous Line 926  class StokesProblemCartesian(Homogeneous
926           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
927           :rtype: ``float``           :rtype: ``float``
928           """           """
929           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
930    
931       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
932           """           """
# Line 937  class StokesProblemCartesian(Homogeneous Line 937  class StokesProblemCartesian(Homogeneous
937           :return: inner product of p0 and p1           :return: inner product of p0 and p1
938           :rtype: ``float``           :rtype: ``float``
939           """           """
940           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
941           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
942           return util.integrate(s0*s1)           return util.integrate(s0*s1)
943    
944       def norm_v(self,v):       def norm_v(self,v):
# Line 978  class StokesProblemCartesian(Homogeneous Line 978  class StokesProblemCartesian(Homogeneous
978          :rtype: equal to the type of p          :rtype: equal to the type of p
979          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
980          """          """
981          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
982    
983       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
984           """           """
# Line 988  class StokesProblemCartesian(Homogeneous Line 988  class StokesProblemCartesian(Homogeneous
988           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
989           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
990           """           """
991           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
992           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
993           return  out           return  out
994    

Legend:
Removed from v.3360  
changed lines
  Added in v.3432

  ViewVC Help
Powered by ViewVC 1.1.26