/[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 3074 by gross, Tue Jul 27 01:47:45 2010 UTC revision 3569 by gross, Thu Sep 1 02:42:36 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 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     :note: The problem is solved in a least squares formulation.     :cvar EVAL: direct pressure gradient evaluation for flux
50       :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                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53     """     """
54         EVAL="EVAL"
55     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="EVAL", 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
62        :type domain: `Domain`        :type domain: `Domain`
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 adaptSubTolerance: switches on automatic subtolerance selection        :param solver: solver method
66        :type adaptSubTolerance: ``bool``        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
67        :param solveForFlux: if True the solver solves for the flux (do not use!)        :param verbose: if ``True`` some information on the iteration progress are printed.
68        :type solveForFlux: ``bool``          :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70          :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.solveForFlux=solveForFlux        self.solver=solver
78        self.useReduced=useReduced        self.useReduced=useReduced
79        self.__adaptSubTolerance=adaptSubTolerance        self.verbose=verbose
80        self.verbose=False        self.l=None
81                self.w=None
82        self.__pde_k=LinearPDESystem(domain)      
       self.__pde_k.setSymmetryOn()  
       if self.useReduced: self.__pde_k.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        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting        if self.solver  == self.EVAL:
88        # self.__pde_l.setSymmetryOn()           self.__pde_v=None
89        # if self.useReduced: self.__pde_l.setReducedOrderOn()       if self.verbose: print "DarcyFlow: simple solver is used."
90        self.setTolerance()  
91        self.setAbsoluteTolerance()        elif self.solver  == self.POST:
92        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))       if util.inf(w)<0.:
93        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))          raise ValueError,"Weighting factor must be non-negative."
94               if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95             self.__pde_v=LinearPDESystem(domain)
96             self.__pde_v.setSymmetryOn()
97             if self.useReduced: self.__pde_v.setReducedOrderOn()
98         self.w=w
99             self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100    
101          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"))
109          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112        
113            
114       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
115          """
116          assigns values to model parameters
117    
118          :param f: volumetic sources/sinks
119          :type f: scalar value on the domain (e.g. `escript.Data`)
120          :param g: flux sources/sinks
121          :type g: vector values on the domain (e.g. `escript.Data`)
122          :param location_of_fixed_pressure: mask for locations where pressure is fixed
123          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
124          :param location_of_fixed_flux:  mask for locations where flux is fixed.
125          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
126          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
127          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
128    
129          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
130          :note: at any point on the boundary of the domain the pressure
131                 (``location_of_fixed_pressure`` >0) or the normal component of the
132                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
133                 is along the *x_i* axis.
134    
135          """
136          if location_of_fixed_pressure!=None:
137               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
138               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
139          if location_of_fixed_flux!=None:
140              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
141              if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
142                
143          if permeability!=None:
144        
145         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
146            
147         if perm.getRank()==0:
148    
149            perm_inv=(1./perm)
150            perm_inv=perm_inv*util.kronecker(self.domain.getDim())
151            perm=perm*util.kronecker(self.domain.getDim())
152            
153            
154         elif perm.getRank()==2:
155            perm_inv=util.inverse(perm)
156         else:
157            raise ValueError,"illegal rank of permeability."
158            
159         self.__permeability=perm
160         self.__permeability_inv=perm_inv
161        
162             #====================
163         self.__pde_p.setValue(A=self.__permeability)
164             if self.solver  == self.EVAL:
165                  pass # no extra work required
166             elif self.solver  == self.POST:
167            k=util.kronecker(self.domain.getDim())
168            self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
169            self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
170             elif self.solver  == self.SMOOTH:
171            self.__pde_v.setValue(D=self.__permeability_inv)
172    
173          if g != None:
174        g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
175        if g.isEmpty():
176              g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
177        else:
178            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
179        self.__g=g
180          if f !=None:
181         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182         if f.isEmpty():      
183              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184         else:
185             if f.getRank()>0: raise ValueError,"illegal rank of f."
186         self.__f=f
187    
188     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
189        """        """
190        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
191        :return: `SolverOptions`        :return: `SolverOptions`
192        """        """
193        return self.__pde_k.getSolverOptions()        if self.__pde_v == None:
194              return None
195          else:
196              return self.__pde_v.getSolverOptions()
197                
198     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
199        """        """
200        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
201        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
202        :param options: `SolverOptions`        :param options: `SolverOptions`
       :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
203        """        """
204        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
205              self.__pde_v.setSolverOptions(options)
206            
207     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
208        """        """
209        Returns the solver options used to solve the pressure problems        Returns the solver options used to solve the pressure problems
         
       *(Q^* K Q)p=-Q^*G*  
         
210        :return: `SolverOptions`        :return: `SolverOptions`
211        """        """
212        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 119  class DarcyFlow(object): Line 214  class DarcyFlow(object):
214     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
215        """        """
216        Sets the solver options used to solve the pressure problems        Sets the solver options used to solve the pressure problems
         
       *(Q^* K Q)p=-Q^*G*  
         
217        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
218                
219        :param options: `SolverOptions`        :param options: `SolverOptions`
220        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221        """        """
222        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
   
    def getSolverOptionsWeighting(self):  
       """  
       Returns the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       :return: `SolverOptions`  
       """  
       return self.__pde_l.getSolverOptions()  
   
    def setSolverOptionsWeighting(self, options=None):  
       """  
       Sets the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       If ``options`` is not present, the options are reset to default  
   
       :param options: `SolverOptions`  
       :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
       """  
       return self.__pde_l.setSolverOptions(options)  
   
   
    def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):  
       """  
       assigns values to model parameters  
         
       :param f: volumetic sources/sinks  
       :type f: scalar value on the domain (e.g. `Data`)  
       :param g: flux sources/sinks  
       :type g: vector values on the domain (e.g. `Data`)  
       :param location_of_fixed_pressure: mask for locations where pressure is fixed  
       :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)  
       :param location_of_fixed_flux:  mask for locations where flux is fixed.  
       :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)  
       :param permeability: permeability tensor. If scalar ``s`` is given the tensor with  
       ``s`` on the main diagonal is used.  
       :type permeability: scalar or tensor values on the domain (e.g. `Data`)  
       :note: the values of parameters which are not set by calling ``setValue`` are not altered.  
       :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)  
       or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal  
       is along the *x_i* axis.  
       """  
       if f !=None:  
      f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))  
      if f.isEmpty():  
         f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))  
      else:  
         if f.getRank()>0: raise ValueError,"illegal rank of f."  
         self.__f=f  
       if g !=None:  
      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      if g.isEmpty():  
         g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
            raise ValueError,"illegal shape of g"  
         self.__g=g  
       if location_of_fixed_pressure!=None:  
            self.__pde_p.setValue(q=location_of_fixed_pressure)  
            #self.__pde_l.setValue(q=location_of_fixed_pressure)  
       if location_of_fixed_flux!=None:  
            self.__pde_k.setValue(q=location_of_fixed_flux)  
               
       if permeability!=None:  
      perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
      if perm.getRank()==0:  
         perm_inv=(1./perm)*util.kronecker(self.domain.getDim())  
         perm=perm*util.kronecker(self.domain.getDim())  
      elif perm.getRank()==2:  
         perm_inv=util.inverse(perm)  
      else:  
         raise ValueError,"illegal rank of permeability."  
   
      self.__permeability=perm  
      self.__permeability_inv=perm_inv  
      self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10  
      #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10  
      if  self.solveForFlux:  
         self.__pde_k.setValue(D=self.__permeability_inv)  
      else:  
         self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
      self.__pde_p.setValue(A=self.__permeability)  
      #self.__pde_l.setValue(D=1/self.__l)  
          #self.__pde_l.setValue(A=self.__permeability)  
   
    def __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l  
       else:  
      return (f-util.div(v))*self.__l  
       # if f == None:  
       #      self.__pde_l.setValue(Y=-util.div(v))    
       # else:  
       #      return (f-util.div(v))/self.__l  
       # return self.__pde_l.getSolution()  
         
    def __getPressure(self, v, p0, g=None):  
       # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0  
       if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"  
       if g == None:  
      self.__pde_p.setValue(X=-v, r=p0)  
       else:  
      self.__pde_p.setValue(X=g-v, r=p0)        
       p=self.__pde_p.getSolution()  
       return p  
   
    def __Aprod_v(self,dv):  
       # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv    
       dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)  
       a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp  
       b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)  
       return ArithmeticTuple(a,b,-dp)  
   
    def __Msolve_PCG_v(self,r):  
       # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv  
       if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"  
       self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_k.getSolution()  
   
    def __inner_PCG_v(self,v,r):  
       return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])  
223                
224     def __Aprod_p(self,dp):     def solve(self, u0, p0):
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())  
       du=self.__pde_k.getSolution()  
       # self.__pde_v.getOperator().saveMM("proj.mm")  
       return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)  
   
    def __getFlux(self,p, v0, f=None, g=None):  
       # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp  
       if f!=None:  
      self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))  
       self.__pde_k.setValue(r=v0)  
       g2=util.tensor_mult(self.__permeability_inv,g)  
       if p == None:  
      self.__pde_k.setValue(Y=g2)  
       else:  
      self.__pde_k.setValue(Y=g2-util.grad(p))  
       return self.__pde_k.getSolution()    
         
       #v=self.__getFlux(p, u0, f=self.__f, g=g2)        
    def __Msolve_PCG_p(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
       self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
   
    def __L2(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
    def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
225        """        """
226        solves the problem.        solves the problem.
227                
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
228        :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.
229        :type u0: vector value on the domain (e.g. `Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
230        :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.
231        :type p0: scalar value on the domain (e.g. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
232        :return: flux and pressure        :return: flux and pressure
233        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.Data`.
   
       :note: The problem is solved as a least squares form  
       *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*  
       *G^*u+G^* K Gp=G^*g*  
   
       where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.  
       """  
       self.verbose=verbose  
       rtol=self.getTolerance()  
       atol=self.getAbsoluteTolerance()  
       self.setSubProblemTolerance()  
       num_corrections=0  
       converged=False  
       norm_r=None  
         
       # Eliminate the hydrostatic pressure:  
       if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."  
       self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))          
       p0=self.__pde_p.getSolution()  
       g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))  
       norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5      
   
       p=p0*0  
       if self.solveForFlux:  
      v=u0.copy()  
       else:  
      v=self.__getFlux(p, u0, f=self.__f, g=g2)  
   
       while not converged and norm_g2 > 0:  
      Gp=util.grad(p)  
      KGp=util.tensor_mult(self.__permeability,Gp)  
      if self.verbose:  
         def_p=g2-(v+KGp)  
         def_v=self.__f-util.div(v)  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))  
         print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5  
         print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2  
         print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5  
      ATOL=atol+rtol*norm_g2  
      if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)  
      if norm_r == None or norm_r>ATOL:  
         if num_corrections>max_num_corrections:  
            raise ValueError,"maximum number of correction steps reached."  
         
         if self.solveForFlux:  
            # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)  
            v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),  
                self.__Aprod_v,  
                v,  
                self.__Msolve_PCG_v,  
                self.__inner_PCG_v,  
                atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            p=r[2]  
         else:  
            # initial residual is r=G^*(g2-KGp - v)  
            p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),  
                  self.__Aprod_p,  
                  p,  
                  self.__Msolve_PCG_p,  
                  self.__inner_PCG_p,  
                  atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            v=r[1]  
         if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r  
         num_corrections+=1  
      else:  
         if self.verbose: print "DarcyFlux: stopping criterium reached."  
         converged=True  
       return v,p+p0  
    def setTolerance(self,rtol=1e-4):  
       """  
       sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
234    
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
         
       where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).  
         
       :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 setAbsoluteTolerance(self,atol=0.):  
235        """        """
236        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        self.__pde_p.setValue(X=self.__g ,
237                                Y=self.__f,
238                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
239                                r=p0)
240          p=self.__pde_p.getSolution()
241          u = self.getFlux(p, u0)
242          return u,p
243                
244        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*     def getFlux(self,p, u0=None):
   
   
       where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
   
       :param atol: absolute tolerance for the pressure  
       :type atol: non-negative ``float``  
       """  
       if atol<0:  
      raise ValueError,"Absolute tolerance needs to be non-negative."  
       self.__atol=atol  
    def getAbsoluteTolerance(self):  
       """  
       returns the absolute tolerance  
       :return: current absolute tolerance  
       :rtype: ``float``  
       """  
       return self.__atol  
    def getSubProblemTolerance(self):  
       """  
       Returns a suitable subtolerance  
       :type: ``float``  
       """  
       return max(util.EPSILON**(0.5),self.getTolerance()**2)  
   
    def setSubProblemTolerance(self):  
       """  
       Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
       """  
       if self.__adaptSubTolerance:  
      sub_tol=self.getSubProblemTolerance()  
      self.getSolverOptionsFlux().setTolerance(sub_tol)  
      self.getSolverOptionsFlux().setAbsoluteTolerance(0.)  
      self.getSolverOptionsPressure().setTolerance(sub_tol)  
      self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
      self.getSolverOptionsWeighting().setTolerance(sub_tol)  
      self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)  
      if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
   
 class DarcyFlowOld(object):  
     """  
     solves the problem  
   
     *u_i+k_{ij}*p_{,j} = g_i*  
     *u_{i,i} = f*  
   
     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,  
   
     :note: The problem is solved in a least squares formulation.  
     """  
   
     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):  
         """  
         initializes the Darcy flux problem  
         :param domain: domain of the problem  
         :type domain: `Domain`  
     :param useReduced: uses reduced oreder on flux and pressure  
     :type useReduced: ``bool``  
     :param adaptSubTolerance: switches on automatic subtolerance selection  
     :type adaptSubTolerance: ``bool``    
245          """          """
246          self.domain=domain          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
         if weight == None:  
            s=self.domain.getSize()  
            self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2  
            # self.__l=(3.*util.longestEdge(self.domain))**2  
            #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2  
         else:  
            self.__l=weight  
         self.__pde_v=LinearPDESystem(domain)  
         if useReduced: self.__pde_v.setReducedOrderOn()  
         self.__pde_v.setSymmetryOn()  
         self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))  
         self.__pde_p=LinearSinglePDE(domain)  
         self.__pde_p.setSymmetryOn()  
         if useReduced: self.__pde_p.setReducedOrderOn()  
         self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
         self.setTolerance()  
         self.setAbsoluteTolerance()  
     self.__adaptSubTolerance=adaptSubTolerance  
     self.verbose=False  
     def getSolverOptionsFlux(self):  
     """  
     Returns the solver options used to solve the flux problems  
       
     *(I+D^*D)u=F*  
       
     :return: `SolverOptions`  
     """  
     return self.__pde_v.getSolverOptions()  
     def setSolverOptionsFlux(self, options=None):  
     """  
     Sets the solver options used to solve the flux problems  
       
     *(I+D^*D)u=F*  
       
     If ``options`` is not present, the options are reset to default  
     :param options: `SolverOptions`  
     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
     """  
     return self.__pde_v.setSolverOptions(options)  
     def getSolverOptionsPressure(self):  
     """  
     Returns the solver options used to solve the pressure problems  
       
     *(Q^*Q)p=Q^*G*  
       
     :return: `SolverOptions`  
     """  
     return self.__pde_p.getSolverOptions()  
     def setSolverOptionsPressure(self, options=None):  
     """  
     Sets the solver options used to solve the pressure problems  
       
     *(Q^*Q)p=Q^*G*  
       
     If ``options`` is not present, the options are reset to default  
     :param options: `SolverOptions`  
     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
     """  
     return self.__pde_p.setSolverOptions(options)  
   
     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):  
         """  
         assigns values to model parameters  
   
         :param f: volumetic sources/sinks  
         :type f: scalar value on the domain (e.g. `Data`)  
         :param g: flux sources/sinks  
         :type g: vector values on the domain (e.g. `Data`)  
         :param location_of_fixed_pressure: mask for locations where pressure is fixed  
         :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)  
         :param location_of_fixed_flux:  mask for locations where flux is fixed.  
         :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)  
         :param permeability: permeability tensor. If scalar ``s`` is given the tensor with  
                              ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with  
                              ``v`` on the main diagonal is used.  
         :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)  
   
         :note: the values of parameters which are not set by calling ``setValue`` are not altered.  
         :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)  
                or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal  
                is along the *x_i* axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.__f=f  
         if g !=None:  
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
   
     def setTolerance(self,rtol=1e-4):  
         """  
         sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
   
         :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 setAbsoluteTolerance(self,atol=0.):  
         """  
         sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
   
         :param atol: absolute tolerance for the pressure  
         :type atol: non-negative ``float``  
         """  
         if atol<0:  
             raise ValueError,"Absolute tolerance needs to be non-negative."  
         self.__atol=atol  
     def getAbsoluteTolerance(self):  
        """  
        returns the absolute tolerance  
         
        :return: current absolute tolerance  
        :rtype: ``float``  
        """  
        return self.__atol  
     def getSubProblemTolerance(self):  
     """  
     Returns a suitable subtolerance  
     @type: ``float``  
     """  
     return max(util.EPSILON**(0.75),self.getTolerance()**2)  
     def setSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
          """  
      if self.__adaptSubTolerance:  
          sub_tol=self.getSubProblemTolerance()  
              self.getSolverOptionsFlux().setTolerance(sub_tol)  
          self.getSolverOptionsFlux().setAbsoluteTolerance(0.)  
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
          """  
          solves the problem.  
   
          The iteration is terminated if the residual norm is less then self.getTolerance().  
   
          :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.  
          :type u0: vector value on the domain (e.g. `Data`).  
          :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.  
          :type p0: scalar value on the domain (e.g. `Data`).  
          :param verbose: if set some information on iteration progress are printed  
          :type verbose: ``bool``  
          :return: flux and pressure  
          :rtype: ``tuple`` of `Data`.  
   
          :note: The problem is solved as a least squares form  
   
          *(I+D^*D)u+Qp=D^*f+g*  
          *Q^*u+Q^*Qp=Q^*g*  
   
          where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
          We eliminate the flux form the problem by setting  
   
          *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux  
   
          form the first equation. Inserted into the second equation we get  
   
          *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure  
   
          which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step  
          PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
      self.setSubProblemTolerance()  
          num_corrections=0  
          converged=False  
          p=p0  
          norm_r=None  
          while not converged:  
                v=self.getFlux(p, fixed_flux=u0)  
                Qp=self.__Q(p)  
                norm_v=self.__L2(v)  
                norm_Qp=self.__L2(Qp)  
                if norm_v == 0.:  
                   if norm_Qp == 0.:  
                      return v,p  
                   else:  
                     fac=norm_Qp  
                else:  
                   if norm_Qp == 0.:  
                     fac=norm_v  
                   else:  
                     fac=2./(1./norm_v+1./norm_Qp)  
                ATOL=(atol+rtol*fac)  
                if self.verbose:  
                     print "DarcyFlux: L2 norm of v = %e."%norm_v  
                     print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp  
                     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 div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)  
                     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
                if norm_r == None or norm_r>ATOL:  
                    if num_corrections>max_num_corrections:  
                          raise ValueError,"maximum number of correction steps reached."  
                    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)  
                    num_corrections+=1  
                else:  
                    converged=True  
          return v,p  
     def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
     def __Q(self,p):  
           return util.tensor_mult(self.__permeability,util.grad(p))  
   
     def __Aprod(self,dp):  
           if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
           Qdp=self.__Q(dp)  
           self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution()  
           # self.__pde_v.getOperator().saveMM("proj.mm")  
           return Qdp+du  
     def __inner_GMRES(self,r,s):  
          return util.integrate(util.inner(r,s))  
   
     def __inner_PCG(self,p,r):  
          return util.integrate(util.inner(self.__Q(p), r))  
   
     def __Msolve_PCG(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())  
           # self.__pde_p.getOperator().saveMM("prec.mm")  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=Data()):  
         """  
         returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``  
247          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
248          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
249    
250          :param p: pressure.          :param p: pressure.
251          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
252          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
253          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
254          :return: flux          :return: flux
255          :rtype: `Data`          :rtype: `escript.Data`
         :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*  
                for the permeability *k_{ij}*  
256          """          """
257      self.setSubProblemTolerance()          u_eval=self.__g-util.tensor_mult(self.__permeability,util.grad(p))
258          g=self.__g          if self.solver  == self.EVAL:
259          f=self.__f             u = self.__g-util.tensor_mult(self.__permeability,util.grad(p))
260          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
261          if p == None:              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p))
262             self.__pde_v.setValue(Y=g)              if u0 == None:
263          else:             self.__pde_v.setValue(r=escript.Data())
264             self.__pde_v.setValue(Y=g-self.__Q(p))          else:
265          return self.__pde_v.getSolution()             self.__pde_v.setValue(r=u0)
266                u= self.__pde_v.getSolution()
267        return u
268          
269  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
270       """       """
271       solves       solves
# Line 792  class StokesProblemCartesian(Homogeneous Line 298  class StokesProblemCartesian(Homogeneous
298           """           """
299           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
300           self.domain=domain           self.domain=domain
301           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
302           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
303            
304           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
305           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 810  class StokesProblemCartesian(Homogeneous Line 316  class StokesProblemCartesian(Homogeneous
316            
317       :rtype: `SolverOptions`       :rtype: `SolverOptions`
318       """       """
319       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
320       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
321           """           """
322       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 818  class StokesProblemCartesian(Homogeneous Line 324  class StokesProblemCartesian(Homogeneous
324       :param options: new solver  options       :param options: new solver  options
325       :type options: `SolverOptions`       :type options: `SolverOptions`
326       """       """
327           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
328       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
329           """           """
330       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 875  class StokesProblemCartesian(Homogeneous Line 381  class StokesProblemCartesian(Homogeneous
381          if eta !=None:          if eta !=None:
382              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
383              kk=util.outer(k,k)              kk=util.outer(k,k)
384              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
385          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
386              self.__pde_u.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)))
387          if restoration_factor!=None:          if restoration_factor!=None:
388              n=self.domain.getNormal()              n=self.domain.getNormal()
389              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
390          if fixed_u_mask!=None:          if fixed_u_mask!=None:
391              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
392          if f!=None: self.__f=f          if f!=None: self.__f=f
393          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
394          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
395    
396       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):
397          """          """
398          assigns values to the model parameters          assigns values to the model parameters
399    
# Line 927  class StokesProblemCartesian(Homogeneous Line 433  class StokesProblemCartesian(Homogeneous
433           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
434           :rtype: ``float``           :rtype: ``float``
435           """           """
436           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)))
437    
438       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
439           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 444  class StokesProblemCartesian(Homogeneous
444           :return: inner product of p0 and p1           :return: inner product of p0 and p1
445           :rtype: ``float``           :rtype: ``float``
446           """           """
447           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
448           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
449           return util.integrate(s0*s1)           return util.integrate(s0*s1)
450    
451       def norm_v(self,v):       def norm_v(self,v):
# Line 962  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
468           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
469           """           """
470           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
471           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
472       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
473       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
474           if self.__stress.isEmpty():           if self.__stress.isEmpty():
475              self.__pde_u.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)))
476           else:           else:
477              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
478           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
479           return  out           return  out
480    
481       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 979  class StokesProblemCartesian(Homogeneous Line 485  class StokesProblemCartesian(Homogeneous
485          :rtype: equal to the type of p          :rtype: equal to the type of p
486          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
487          """          """
488          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))
489    
490       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
491           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 495  class StokesProblemCartesian(Homogeneous
495           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
496           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
497           """           """
498           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
499           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
500           return  out           return  out
501    
502       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

Legend:
Removed from v.3074  
changed lines
  Added in v.3569

  ViewVC Help
Powered by ViewVC 1.1.26