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

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

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

revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC revision 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC
# Line 1  Line 1 
1  # -*- coding: utf-8 -*-  # -*- coding: utf-8 -*-
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2012 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # http://www.uq.edu.au  # http://www.uq.edu.au
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
# Line 13  Line 13 
13  #  #
14  ##############################################################################  ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2012 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17  http://www.uq.edu.au  http://www.uq.edu.au
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
# Line 33  Some models for flow Line 33  Some models for flow
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36  from . import escript  from . import escriptcpp
37    escore=escriptcpp
38    #from . import escript
39  from . import util  from . import util
40  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
41  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
# Line 58  class DarcyFlow(object): Line 60  class DarcyFlow(object):
60     SMOOTH="SMOOTH"     SMOOTH="SMOOTH"
61     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
62        """        """
63        initializes the Darcy flux problem        initializes the Darcy flux problem.
64    
65        :param domain: domain of the problem        :param domain: domain of the problem
66        :type domain: `Domain`        :type domain: `Domain`
67        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
68        :type useReduced: ``bool``        :type useReduced: ``bool``
69        :param solver: solver method        :param solver: solver method
70        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
71        :param verbose: if ``True`` some information on the iteration progress are printed.        :param verbose: if ``True`` some information on the iteration progress are printed.
72        :type verbose: ``bool``        :type verbose: ``bool``
73        :param w: weighting factor for `DarcyFlow.POST` solver        :param w: weighting factor for `DarcyFlow.POST` solver
# Line 98  class DarcyFlow(object): Line 101  class DarcyFlow(object):
101           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
102           self.w=w           self.w=w
103           x=self.domain.getX()           x=self.domain.getX()
104           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in range(self.domain.getDim()) ] )
105           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
106    
107        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
# Line 108  class DarcyFlow(object): Line 111  class DarcyFlow(object):
111           if self.verbose: print("DarcyFlow: flux smoothing is used.")           if self.verbose: print("DarcyFlow: flux smoothing is used.")
112           self.w=0           self.w=0
113    
114        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
115        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
116        self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
117        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118        self.ref_point_id=None        self.ref_point_id=None
119        self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)        self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
120        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
121        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_flux = escore.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
122        self.perm_scale=1.        self.perm_scale=1.
123            
124                    
# Line 288  class DarcyFlow(object): Line 291  class DarcyFlow(object):
291          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
292          return self._getFlux(p-p_hydrostatic, u0)          return self._getFlux(p-p_hydrostatic, u0)
293    
294     def _getFlux(self,pp, u0=None):     def _getFlux(self, pp, u0=None):
295          """          """
296          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``          returns the flux for a given pressure ``pp`` where the flux is equal to
297          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
298          Notice that ``g`` is used, see `setValue`.          `setValue`). Notice that ``g`` is used, see `setValue`.
299    
300          :param p: pressure.          :param pp: pressure.
301          :type p: scalar value on the domain (e.g. `escript.Data`).          :type pp: scalar value on the domain (i.e. `escript.Data`).
302          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
303          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``          :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
304          :return: flux          :return: flux
305          :rtype: `escript.Data`          :rtype: `escript.Data`
306          """          """
# Line 307  class DarcyFlow(object): Line 310  class DarcyFlow(object):
310              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
311              print              print
312              if u0 == None:              if u0 == None:
313                 self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escore.Data())
314              else:              else:
315                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))                 if not isinstance(u0, escore.Data) : u0 = escore.Vector(u0, escore.Solution(self.domain))
316                 self.__pde_v.setValue(r=1./self.perm_scale * u0)                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
317              u= self.__pde_v.getSolution() * self.perm_scale              u= self.__pde_v.getSolution() * self.perm_scale
318          return u          return u
# Line 431  class StokesProblemCartesian(Homogeneous Line 434  class StokesProblemCartesian(Homogeneous
434          if eta !=None:          if eta !=None:
435              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
436              kk=util.outer(k,k)              kk=util.outer(k,k)
437              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escore.Function(self.domain))
438              self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
439              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
440          if restoration_factor!=None:          if restoration_factor!=None:
# Line 443  class StokesProblemCartesian(Homogeneous Line 446  class StokesProblemCartesian(Homogeneous
446          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
447          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
448    
449       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):       def initialize(self,f=escore.Data(),fixed_u_mask=escore.Data(),eta=1,surface_stress=escore.Data(),stress=escore.Data(), restoration_factor=0):
450          """          """
451          assigns values to the model parameters          assigns values to the model parameters
452    
# Line 483  class StokesProblemCartesian(Homogeneous Line 486  class StokesProblemCartesian(Homogeneous
486           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
487           :rtype: ``float``           :rtype: ``float``
488           """           """
489           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))           return util.integrate(util.interpolate(p,escore.Function(self.domain))*util.interpolate(Bv, escore.Function(self.domain)))
490    
491       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
492           """           """
# Line 494  class StokesProblemCartesian(Homogeneous Line 497  class StokesProblemCartesian(Homogeneous
497           :return: inner product of p0 and p1           :return: inner product of p0 and p1
498           :rtype: ``float``           :rtype: ``float``
499           """           """
500           s0=util.interpolate(p0, escript.Function(self.domain))           s0=util.interpolate(p0, escore.Function(self.domain))
501           s1=util.interpolate(p1, escript.Function(self.domain))           s1=util.interpolate(p1, escore.Function(self.domain))
502           return util.integrate(s0*s1)           return util.integrate(s0*s1)
503    
504       def norm_v(self,v):       def norm_v(self,v):
# Line 535  class StokesProblemCartesian(Homogeneous Line 538  class StokesProblemCartesian(Homogeneous
538          :rtype: equal to the type of p          :rtype: equal to the type of p
539          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
540          """          """
541          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escore.Function(self.domain))**2))
542    
543       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
544           """           """
# Line 545  class StokesProblemCartesian(Homogeneous Line 548  class StokesProblemCartesian(Homogeneous
548           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
549           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
550           """           """
551           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escore.Data(), y=escore.Data(), X=-p*util.kronecker(self.domain))
552           out=self.__pde_v.getSolution()           out=self.__pde_v.getSolution()
553           return  out           return  out
554    
555       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
556           """           """
557           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
558           with accuracy `self.getSubProblemTolerance()`           with accuracy ``self.getSubProblemTolerance()``
559    
560           :param Bv: velocity increment           :param Bv: velocity increment
561           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*

Legend:
Removed from v.3981  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26