/[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 3852 by jfenwick, Thu Mar 1 05:34:52 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-2010 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2010 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
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
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 32  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 57  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 96  class DarcyFlow(object): Line 100  class DarcyFlow(object):
100           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
101           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
102           self.w=w           self.w=w
103           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           x=self.domain.getX()
104             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
106    
107        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
108           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
# Line 105  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.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
117        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118          self.ref_point_id=None
119          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
120          self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
121          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 135  class DarcyFlow(object): Line 145  class DarcyFlow(object):
145    
146        """        """
147        if location_of_fixed_pressure!=None:        if location_of_fixed_pressure!=None:
148             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
149               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
150               if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
151               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
152               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
153             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
154        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
155            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
156            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)            if not self.__pde_v == None:
157                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
158                            
159        if permeability!=None:        if permeability!=None:
160            
# Line 170  class DarcyFlow(object): Line 185  class DarcyFlow(object):
185           elif self.solver  == self.POST:           elif self.solver  == self.POST:
186                k=util.kronecker(self.domain.getDim())                k=util.kronecker(self.domain.getDim())
187                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
188                self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
189                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
190           elif self.solver  == self.SMOOTH:           elif self.solver  == self.SMOOTH:
191              self.__pde_v.setValue(D=self.__permeability_inv)              self.__pde_v.setValue(D=self.__permeability_inv)
192    
# Line 181  class DarcyFlow(object): Line 197  class DarcyFlow(object):
197          else:          else:
198               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
199          self.__g=g          self.__g=g
200            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
201            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
202        if f !=None:        if f !=None:
203           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204           if f.isEmpty():                 if f.isEmpty():      
# Line 237  class DarcyFlow(object): Line 255  class DarcyFlow(object):
255        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
256    
257        """        """
258        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
259          if self.ref_point_id == None:
260              p_ref=0
261          else:
262              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
263          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
264          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
265          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
266                              Y=self.__f * 1./self.perm_scale,                              Y=self.__f * 1./self.perm_scale,
267                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
268                              r=p0)                              r=p0 - p0_hydrostatic)
269        p=self.__pde_p.getSolution()        pp=self.__pde_p.getSolution()
270        u = self.getFlux(p, u0)        u = self._getFlux(pp, u0)
271        return u,p        return u,pp + p0_hydrostatic
272                
273     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
274          """          """
# Line 258  class DarcyFlow(object): Line 283  class DarcyFlow(object):
283          :return: flux          :return: flux
284          :rtype: `escript.Data`          :rtype: `escript.Data`
285          """          """
286            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
287            if self.ref_point_id == None:
288                p_ref=0
289            else:
290                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
291            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)
293    
294       def _getFlux(self, pp, u0=None):
295            """
296            returns the flux for a given pressure ``pp`` where the flux is equal to
297            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
298            `setValue`). Notice that ``g`` is used, see `setValue`.
299    
300            :param pp: pressure.
301            :type pp: scalar value on the domain (i.e. `escript.Data`).
302            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
303            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
304            :return: flux
305            :rtype: `escript.Data`
306            """
307          if self.solver  == self.EVAL:          if self.solver  == self.EVAL:
308             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
309          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
310              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
311                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
319                
320  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
# Line 387  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 399  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 439  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 450  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 491  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 501  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.3852  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26