/[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 3515 by gross, Thu May 19 08:20:57 2011 UTC revision 3885 by gross, Wed Apr 4 22:12:26 2012 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  import escript  from . import escript
36  import util  from . 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
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41     """     """
# Line 46  class DarcyFlow(object): Line 46  class DarcyFlow(object):
46        
47     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48        
49     :cvar SIMPLE: simple solver     :cvar EVAL: direct pressure gradient evaluation for flux
50     :cvar POST: solver using global postprocessing of flux     :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
51     :cvar STAB: solver uses (non-symmetric) stabilization                 where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52     :cvar SYMSTAB: solver uses symmetric stabilization     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53     """     """
54     SIMPLE="SIMPLE"     EVAL="EVAL"
55       SIMPLE="EVAL"
56     POST="POST"     POST="POST"
57     STAB="STAB"     SMOOTH="SMOOTH"
58     SYMSTAB="SYMSTAB"     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
    def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):  
59        """        """
60        initializes the Darcy flux problem        initializes the Darcy flux problem
61        :param domain: domain of the problem        :param domain: domain of the problem
# Line 63  class DarcyFlow(object): Line 63  class DarcyFlow(object):
63        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
64        :type useReduced: ``bool``        :type useReduced: ``bool``
65        :param solver: solver method        :param solver: solver method
66        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
67        :param verbose: if ``True`` some information on the iteration progress are printed.        :param verbose: if ``True`` some information on the iteration progress are printed.
68        :type verbose: ``bool``        :type verbose: ``bool``
69        :param w: weighting factor for `DarcyFlow.POST` solver        :param w: weighting factor for `DarcyFlow.POST` solver
70        :type w: ``float``        :type w: ``float``
71                
72        """        """
73          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError("unknown solver %d."%solver)
75    
76        self.domain=domain        self.domain=domain
77        self.solver=solver        self.solver=solver
78        self.useReduced=useReduced        self.useReduced=useReduced
79        self.verbose=verbose        self.verbose=verbose
80        self.scale=1.        self.l=None
81                self.w=None
82              
       self.__pde_v=LinearPDESystem(domain)  
       self.__pde_v.setSymmetryOn()  
       if self.useReduced: self.__pde_v.setReducedOrderOn()  
83        self.__pde_p=LinearSinglePDE(domain)        self.__pde_p=LinearSinglePDE(domain)
84        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
85        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
86          
87        if self.solver  == self.SIMPLE:        if self.solver  == self.EVAL:
88       if self.verbose: print "DarcyFlow: simple solver is used."           self.__pde_v=None
89           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))           if self.verbose: print("DarcyFlow: simple solver is used.")
90    
91        elif self.solver  == self.POST:        elif self.solver  == self.POST:
92       self.w=w           if util.inf(w)<0.:
93       if util.inf(w)<0.:              raise ValueError("Weighting factor must be non-negative.")
94          raise ValueError,"Weighting factor must be non-negative."           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
95       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."           self.__pde_v=LinearPDESystem(domain)
96        elif self.solver  == self.STAB:           self.__pde_v.setSymmetryOn()
97        if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."           if self.useReduced: self.__pde_v.setReducedOrderOn()
98        elif  self.solver  == self.SYMSTAB:           self.w=w
99        if self.verbose: print "DarcyFlow: symmetric stabilization is used."           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100        else:  
101      raise ValueError,"unknown solver %s"%self.solver        elif self.solver  == self.SMOOTH:
102             self.__pde_v=LinearPDESystem(domain)
103             self.__pde_v.setSymmetryOn()
104             if self.useReduced: self.__pde_v.setReducedOrderOn()
105             if self.verbose: print("DarcyFlow: flux smoothing is used.")
106             self.w=0
107    
108        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111        self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112        self.setTolerance()        self.perm_scale=1.
113            
114                    
115     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):
# Line 132  class DarcyFlow(object): Line 139  class DarcyFlow(object):
139             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
141            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142            self.__pde_v.setValue(q=self.location_of_fixed_flux)            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
         
143                            
       # pressure  is rescaled by the factor 1/self.scale  
144        if permeability!=None:        if permeability!=None:
145            
146       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147           V=util.vol(self.domain)           self.perm_scale=util.Lsup(util.length(perm))
148           l=V**(1./self.domain.getDim())           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149             perm=perm*(1./self.perm_scale)
150                    
151       if perm.getRank()==0:           if perm.getRank()==0:
152          perm_inv=(1./perm)  
153              self.scale=util.integrate(perm_inv)/V*l              perm_inv=(1./perm)
154          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))              perm=perm*util.kronecker(self.domain.getDim())
156                    
157                    
158       elif perm.getRank()==2:           elif perm.getRank()==2:
159          perm_inv=util.inverse(perm)              perm_inv=util.inverse(perm)
160              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l           else:
161          perm_inv*=(1./self.scale)              raise ValueError("illegal rank of permeability.")
         perm=perm*self.scale  
      else:  
         raise ValueError,"illegal rank of permeability."  
162                    
163       self.__permeability=perm           self.__permeability=perm
164       self.__permeability_inv=perm_inv           self.__permeability_inv=perm_inv
      if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale  
165            
166       if self.solver  == self.SIMPLE:           #====================
167          self.__pde_p.setValue(A=self.__permeability)           self.__pde_p.setValue(A=self.__permeability)
168       elif self.solver  == self.POST:           if self.solver  == self.EVAL:
169          self.__pde_p.setValue(A=self.__permeability)                pass # no extra work required
170          k=util.kronecker(self.domain.getDim())           elif self.solver  == self.POST:
171          self.lamb = self.w*util.length(perm_inv)*l                k=util.kronecker(self.domain.getDim())
172          self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
173       elif self.solver  == self.STAB:                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
174          self.__pde_p.setValue(A=0.5*self.__permeability)                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
175          self.__pde_v.setValue(D=0.5*self.__permeability_inv)           elif self.solver  == self.SMOOTH:
176       elif  self.solver  == self.SYMSTAB:              self.__pde_v.setValue(D=self.__permeability_inv)
         self.__pde_p.setValue(A=0.5*self.__permeability)  
         self.__pde_v.setValue(D=0.5*self.__permeability_inv)  
177    
178        if g != None:        if g != None:
179      g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
180      if g.isEmpty():          if g.isEmpty():
181            g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182      else:          else:
183          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")
184      self.__g=g          self.__g=g
185        if f !=None:        if f !=None:
186       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187       if f.isEmpty():                 if f.isEmpty():      
188            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189       else:           else:
190           if f.getRank()>0: raise ValueError,"illegal rank of f."               if f.getRank()>0: raise ValueError("illegal rank of f.")
191       self.__f=f           self.__f=f
192    
193     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
194        """        """
195        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
196        :return: `SolverOptions`        :return: `SolverOptions`
197        """        """
198        return self.__pde_v.getSolverOptions()        if self.__pde_v == None:
199              return None
200          else:
201              return self.__pde_v.getSolverOptions()
202                
203     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
204        """        """
# Line 202  class DarcyFlow(object): Line 206  class DarcyFlow(object):
206        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
207        :param options: `SolverOptions`        :param options: `SolverOptions`
208        """        """
209        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
210              self.__pde_v.setSolverOptions(options)
211            
212     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
213        """        """
# Line 221  class DarcyFlow(object): Line 226  class DarcyFlow(object):
226        """        """
227        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
228                
229     def setTolerance(self,rtol=1e-4):     def solve(self, u0, p0):
       """  
       sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.  
         
       :param rtol: relative tolerance for the pressure  
       :type rtol: non-negative ``float``  
       """  
       if rtol<0:  
      raise ValueError,"Relative tolerance needs to be non-negative."  
       self.__rtol=rtol  
         
    def getTolerance(self):  
       """  
       returns the relative tolerance  
       :return: current relative tolerance  
       :rtype: ``float``  
       """  
       return self.__rtol  
         
    def solve(self,u0,p0, max_iter=100, iter_restart=20):  
230        """        """
231        solves the problem.        solves the problem.
232                
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
233        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
234        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
235        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
236        :type p0: scalar value on the domain (e.g. `escript.Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
       :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,  
       :type max_iter: ``int``  
       :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.  
                            A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure  
                            of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.  
       :type iter_restart: ``int``  
237        :return: flux and pressure        :return: flux and pressure
238        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
239    
240        """        """
241        # rescale initial guess:        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
242        p0=p0/self.scale                              Y=self.__f * 1./self.perm_scale,
243        if self.solver  == self.SIMPLE or self.solver  == self.POST :                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
244          self.__pde_p.setValue(X=self.__g ,                              r=p0)
245                                Y=self.__f,        p=self.__pde_p.getSolution()
246                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),        u = self.getFlux(p, u0)
                               r=p0)  
         p=self.__pde_p.getSolution()  
         u = self.getFlux(p, u0)  
       elif  self.solver  == self.STAB:  
     u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)  
       elif  self.solver  == self.SYMSTAB:  
     u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)  
       
       if self.verbose:  
         KGp=util.tensor_mult(self.__permeability,util.grad(p))  
         def_p=self.__g-(u+KGp)  
         def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))  
         print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))  
       #rescale result  
       p=p*self.scale  
247        return u,p        return u,p
248                
249     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
250          """          """
251          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
252          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
253          Notice that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
254    
255          :param p: pressure.          :param p: pressure.
256          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
# Line 297  class DarcyFlow(object): Line 259  class DarcyFlow(object):
259          :return: flux          :return: flux
260          :rtype: `escript.Data`          :rtype: `escript.Data`
261          """          """
262          if self.solver  == self.SIMPLE or self.solver  == self.POST  :          if self.solver  == self.EVAL:
263              KGp=util.tensor_mult(self.__permeability,util.grad(p))             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
264              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
265              if u0 == None:              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
            self.__pde_v.setValue(r=escript.Data())  
         else:  
            self.__pde_v.setValue(r=u0)  
             u= self.__pde_v.getSolution()  
     elif self.solver  == self.POST:  
             self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),  
                                   X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))  
266              if u0 == None:              if u0 == None:
267             self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escript.Data())
268          else:              else:
269             self.__pde_v.setValue(r=u0)                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
270              u= self.__pde_v.getSolution()                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
271      elif self.solver  == self.STAB:                 u= self.__pde_v.getSolution() * self.perm_scale
272           gp=util.grad(p)          return u
          self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),  
                                X= p * util.kronecker(self.domain.getDim()),  
                                y= - p * self.domain.getNormal())                            
          if u0 == None:  
            self.__pde_v.setValue(r=escript.Data())  
          else:  
            self.__pde_v.setValue(r=u0)  
          u= self.__pde_v.getSolution()  
     elif  self.solver  == self.SYMSTAB:  
          gp=util.grad(p)  
          self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),  
                                X= escript.Data() ,  
                                y= escript.Data() )                            
          if u0 == None:  
            self.__pde_v.setValue(r=escript.Data())  
          else:  
            self.__pde_v.setValue(r=u0)  
          u= self.__pde_v.getSolution()  
     return u  
         
       
    def __solve_STAB(self, u0, p0, max_iter, iter_restart):  
           # p0 is used as an initial guess  
       u=self.getFlux(p0, u0)    
           self.__pde_p.setValue( Y=self.__f-util.div(u),  
                                  X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),  
                                  y= escript.Data(),  
                                  r=escript.Data())  
   
       dp=self.__pde_p.getSolution()  
       p=GMRES(dp,  
               self.__STAB_Aprod,  
           p0,  
           self.__inner,  
           atol=self.__norm(p0+dp)*self.getTolerance() ,  
           rtol=0.,  
           iter_max=max_iter,  
           iter_restart=iter_restart,  
           verbose=self.verbose,P_R=None)  
               
           u=self.getFlux(p, u0)  
           return u,p  
   
    def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):  
           # p0 is used as an initial guess  
       u=self.getFlux(p0, u0)    
           self.__pde_p.setValue( Y= self.__f,  
                                  X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),  
                                  y=  -  util.inner(self.domain.getNormal(), u),  
                                  r=escript.Data())  
       dp=self.__pde_p.getSolution()  
273                
       p=GMRES(dp,  
               self.__SYMSTAB_Aprod,  
           p0,  
           self.__inner,  
           atol=self.__norm(p0+dp)*self.getTolerance() ,  
           rtol=0.,  
           iter_max=max_iter,  
           iter_restart=iter_restart,  
           verbose=self.verbose,P_R=None)  
               
           u=self.getFlux(p, u0)  
           return u,p  
   
    def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))        
     
    def __norm(self,r):  
          return util.sqrt(self.__inner(r,r))  
           
    def __inner(self,r,s):  
          return util.integrate(util.inner(r,s), escript.Function(self.domain))  
           
    def __STAB_Aprod(self,p):  
       gp=util.grad(p)  
       self.__pde_v.setValue(Y=-0.5*gp,  
                             X=-p*util.kronecker(self.__pde_v.getDomain()),  
                             y= p * self.domain.getNormal(),    
                             r=escript.Data())  
       u = -self.__pde_v.getSolution()  
       self.__pde_p.setValue(Y=util.div(u),  
                             X=0.5*(u+util.tensor_mult(self.__permeability,gp)),  
                             y=escript.Data(),  
                             r=escript.Data())  
       
       return  self.__pde_p.getSolution()  
     
    def __SYMSTAB_Aprod(self,p):  
       gp=util.grad(p)  
       self.__pde_v.setValue(Y=0.5*gp ,  
                             X=escript.Data(),  
                             y=escript.Data(),    
                             r=escript.Data())  
       u = -self.__pde_v.getSolution()  
       self.__pde_p.setValue(Y=escript.Data(),  
                             X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),  
                             y=escript.Data(),  
                             r=escript.Data())  
       
       return  self.__pde_p.getSolution()  
         
   
274  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
275       """       """
276       solves       solves
# Line 436  class StokesProblemCartesian(Homogeneous Line 289  class StokesProblemCartesian(Homogeneous
289              sp.setTolerance()              sp.setTolerance()
290              sp.initialize(...)              sp.initialize(...)
291              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
292                sp.setStokesEquation(...) # new values for some parameters
293                v1,p1=sp.solve(v,p)
294       """       """
295       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
296           """           """
# Line 459  class StokesProblemCartesian(Homogeneous Line 314  class StokesProblemCartesian(Homogeneous
314    
315           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
316           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
317       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
318           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
319    
320       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 468  class StokesProblemCartesian(Homogeneous Line 323  class StokesProblemCartesian(Homogeneous
323            
324       :rtype: `SolverOptions`       :rtype: `SolverOptions`
325       """       """
326       return self.__pde_v.getSolverOptions()           return self.__pde_v.getSolverOptions()
327       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
328           """           """
329       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 482  class StokesProblemCartesian(Homogeneous Line 337  class StokesProblemCartesian(Homogeneous
337       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
338       :rtype: `SolverOptions`       :rtype: `SolverOptions`
339       """       """
340       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
341       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
342           """           """
343       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
344       :param options: new solver  options       :param options: new solver  options
345       :type options: `SolverOptions`       :type options: `SolverOptions`
346       """       """
347       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
348    
349       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
350           """           """
# Line 499  class StokesProblemCartesian(Homogeneous Line 354  class StokesProblemCartesian(Homogeneous
354       :param options: new solver options       :param options: new solver options
355       :type options: `SolverOptions`       :type options: `SolverOptions`
356       """       """
357       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
358       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
359           """           """
360       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 507  class StokesProblemCartesian(Homogeneous Line 362  class StokesProblemCartesian(Homogeneous
362            
363       :rtype: `SolverOptions`       :rtype: `SolverOptions`
364       """       """
365       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
366    
367       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
368           """           """
369           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
370           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
371           """           """
372           pass           pass
373       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
# Line 534  class StokesProblemCartesian(Homogeneous Line 389  class StokesProblemCartesian(Homogeneous
389              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
390              kk=util.outer(k,k)              kk=util.outer(k,k)
391              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
392          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
393              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)))
394          if restoration_factor!=None:          if restoration_factor!=None:
395              n=self.domain.getNormal()              n=self.domain.getNormal()
# Line 558  class StokesProblemCartesian(Homogeneous Line 413  class StokesProblemCartesian(Homogeneous
413          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
414          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
415          :param stress: initial stress          :param stress: initial stress
416      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
417          """          """
418          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
419    
# Line 571  class StokesProblemCartesian(Homogeneous Line 426  class StokesProblemCartesian(Homogeneous
426           :rtype: ``float``           :rtype: ``float``
427           """           """
428           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
429       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
430       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
431           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
432           return out           return out
433    
# Line 613  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
468    
469       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
470           """           """
471           return the value for v for a given p (overwrite)           return the value for v for a given p
472    
473           :param p: a pressure           :param p: a pressure
474           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
# Line 621  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476           """           """
477           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
478           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
479       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
480       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
481           if self.__stress.isEmpty():           if self.__stress.isEmpty():
482              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
483           else:           else:
# Line 661  class StokesProblemCartesian(Homogeneous Line 516  class StokesProblemCartesian(Homogeneous
516           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
517           """           """
518           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
519       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
520       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
521           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
522           return out           return out

Legend:
Removed from v.3515  
changed lines
  Added in v.3885

  ViewVC Help
Powered by ViewVC 1.1.26