/[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 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC revision 4446 by caltinay, Tue Jun 11 04:00:15 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 as escore
37  from . import util  from . import util
38  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from . import linearPDEs as lpe
39  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from . import pdetools as pdt
40    
41  class DarcyFlow(object):  class DarcyFlow(object):
42     """     """
# Line 82  class DarcyFlow(object): Line 82  class DarcyFlow(object):
82        self.l=None        self.l=None
83        self.w=None        self.w=None
84            
85        self.__pde_p=LinearSinglePDE(domain)        self.__pde_p=lpe.LinearSinglePDE(domain)
86        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
87        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
88    
# Line 94  class DarcyFlow(object): Line 94  class DarcyFlow(object):
94           if util.inf(w)<0.:           if util.inf(w)<0.:
95              raise ValueError("Weighting factor must be non-negative.")              raise ValueError("Weighting factor must be non-negative.")
96           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=lpe.LinearPDESystem(domain)
98           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
99           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
100           self.w=w           self.w=w
101           x=self.domain.getX()           x=self.domain.getX()
102           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()) ] )
103           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104    
105        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
106           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=lpe.LinearPDESystem(domain)
107           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
108           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
109           if self.verbose: print("DarcyFlow: flux smoothing is used.")           if self.verbose: print("DarcyFlow: flux smoothing is used.")
110           self.w=0           self.w=0
111    
112        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114        self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
115        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)
116        self.ref_point_id=None        self.ref_point_id=None
117        self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)        self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118        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"))
119        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"))
120        self.perm_scale=1.        self.perm_scale=1.
121            
122                    
# Line 306  class DarcyFlow(object): Line 306  class DarcyFlow(object):
306             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308              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))
309              print  
310              if u0 == None:              if u0 == None:
311                 self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escore.Data())
312              else:              else:
313                 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))
314                 self.__pde_v.setValue(r=1./self.perm_scale * u0)                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
315              u= self.__pde_v.getSolution() * self.perm_scale              u= self.__pde_v.getSolution() * self.perm_scale
316          return u          return u
317                
318  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(pdt.HomogeneousSaddlePointProblem):
319       """       """
320       solves       solves
321    
# Line 347  class StokesProblemCartesian(Homogeneous Line 347  class StokesProblemCartesian(Homogeneous
347           :param domain: domain of the problem.           :param domain: domain of the problem.
348           :type domain: `Domain`           :type domain: `Domain`
349           """           """
350           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           pdt.HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=lpe.LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
353           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
354            
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=lpe.LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
357           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=lpe.LinearPDE(domain)
360           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
361           self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
# Line 432  class StokesProblemCartesian(Homogeneous Line 432  class StokesProblemCartesian(Homogeneous
432          if eta !=None:          if eta !=None:
433              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
434              kk=util.outer(k,k)              kk=util.outer(k,k)
435              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escore.Function(self.domain))
436              self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
437              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)))
438          if restoration_factor!=None:          if restoration_factor!=None:
# Line 444  class StokesProblemCartesian(Homogeneous Line 444  class StokesProblemCartesian(Homogeneous
444          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
445          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
446    
447       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):
448          """          """
449          assigns values to the model parameters          assigns values to the model parameters
450    
# Line 484  class StokesProblemCartesian(Homogeneous Line 484  class StokesProblemCartesian(Homogeneous
484           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
485           :rtype: ``float``           :rtype: ``float``
486           """           """
487           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)))
488    
489       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
490           """           """
# Line 495  class StokesProblemCartesian(Homogeneous Line 495  class StokesProblemCartesian(Homogeneous
495           :return: inner product of p0 and p1           :return: inner product of p0 and p1
496           :rtype: ``float``           :rtype: ``float``
497           """           """
498           s0=util.interpolate(p0, escript.Function(self.domain))           s0=util.interpolate(p0, escore.Function(self.domain))
499           s1=util.interpolate(p1, escript.Function(self.domain))           s1=util.interpolate(p1, escore.Function(self.domain))
500           return util.integrate(s0*s1)           return util.integrate(s0*s1)
501    
502       def norm_v(self,v):       def norm_v(self,v):
# Line 536  class StokesProblemCartesian(Homogeneous Line 536  class StokesProblemCartesian(Homogeneous
536          :rtype: equal to the type of p          :rtype: equal to the type of p
537          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
538          """          """
539          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))
540    
541       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
542           """           """
# Line 546  class StokesProblemCartesian(Homogeneous Line 546  class StokesProblemCartesian(Homogeneous
546           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
547           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
548           """           """
549           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))
550           out=self.__pde_v.getSolution()           out=self.__pde_v.getSolution()
551           return  out           return  out
552    

Legend:
Removed from v.3990  
changed lines
  Added in v.4446

  ViewVC Help
Powered by ViewVC 1.1.26