/[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 3452 by caltinay, Tue Jan 25 01:53:57 2011 UTC revision 3852 by jfenwick, Thu Mar 1 05:34:52 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
   
 print dir(escript)  
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41     """     """
# Line 48  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, useVPIteration=True, weighting_scale=0.1):     SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="POST", 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 useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.        :param w: weighting factor for `DarcyFlow.POST` solver
70        :type useVPIteration: ``bool``            :type w: ``float``
       """  
       self.domain=domain  
       self.useVPIteration=useVPIteration  
       self.useReduced=useReduced  
       self.weighting_scale=weighting_scale  
       if self.useVPIteration:  
          self.solveForFlux=solveForFlux  
          self.__adaptSubTolerance=adaptSubTolerance  
          self.verbose=False  
           
          self.__pde_k=LinearPDESystem(domain)  
          self.__pde_k.setSymmetryOn()  
          if self.useReduced: self.__pde_k.setReducedOrderOn()  
     
          self.__pde_p=LinearSinglePDE(domain)  
          self.__pde_p.setSymmetryOn()  
          if self.useReduced: self.__pde_p.setReducedOrderOn()  
          self.setTolerance()  
          self.setAbsoluteTolerance()  
       else:  
          self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)  
          self.__pde_k.setSymmetryOn()  
          if self.useReduced: self.__pde_k.setReducedOrderOn()  
          C=self.__pde_k.createCoefficient("C")  
          B=self.__pde_k.createCoefficient("B")  
          for i in range(self.domain.getDim()):  
             C[i,self.domain.getDim(),i]=1  
             B[self.domain.getDim(),i,i]=1  
          self.__pde_k.setValue(C=C, B=B)  
       self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))  
       self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))  
         
    def getSolverOptionsFlux(self):  
       """  
       Returns the solver options used to solve the flux problems  
         
       *K^{-1} u=F*  
         
       :return: `SolverOptions`  
       """  
       return self.__pde_k.getSolverOptions()  
         
    def setSolverOptionsFlux(self, options=None):  
       """  
       Sets the solver options used to solve the flux problems  
         
       *K^{-1}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^* K Q)p=-Q^*G*  
         
       :return: `SolverOptions`  
       """  
       return self.__pde_p.getSolverOptions()  
71                
    def setSolverOptionsPressure(self, options=None):  
72        """        """
73        Sets the solver options used to solve the pressure problems        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74                    raise ValueError("unknown solver %d."%solver)
       *(Q^* K 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)  
   
75    
76          self.domain=domain
77          self.solver=solver
78          self.useReduced=useReduced
79          self.verbose=verbose
80          self.l=None
81          self.w=None
82        
83          self.__pde_p=LinearSinglePDE(domain)
84          self.__pde_p.setSymmetryOn()
85          if self.useReduced: self.__pde_p.setReducedOrderOn()
86    
87          if self.solver  == self.EVAL:
88             self.__pde_v=None
89             if self.verbose: print("DarcyFlow: simple solver is used.")
90    
91          elif self.solver  == self.POST:
92             if util.inf(w)<0.:
93                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          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):
116        """        """
117        assigns values to model parameters        assigns values to model parameters
# Line 156  class DarcyFlow(object): Line 125  class DarcyFlow(object):
125        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
126        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
127        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
128        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
129    
130        :note: the values of parameters which are not set by calling ``setValue`` are not altered.        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
131        :note: at any point on the boundary of the domain the pressure        :note: at any point on the boundary of the domain the pressure
# Line 165  class DarcyFlow(object): Line 134  class DarcyFlow(object):
134               is along the *x_i* axis.               is along the *x_i* axis.
135    
136        """        """
137        if self.useVPIteration:        if location_of_fixed_pressure!=None:
138           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
139           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140        else:        if location_of_fixed_flux!=None:
141           q=self.__pde_k.getCoefficient("q")            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
          if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure  
          if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux  
          self.__pde_k.setValue(q=q)  
143                            
       # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)  
144        if permeability!=None:        if permeability!=None:
145       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))      
146           V=util.vol(self.domain)           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147       if perm.getRank()==0:           self.perm_scale=util.Lsup(util.length(perm))
148          perm_inv=(1./perm)           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149              if self.useVPIteration:           perm=perm*(1./self.perm_scale)
150                self.scale=1.          
151              else:           if perm.getRank()==0:
               self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)  
152    
153          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))              perm_inv=(1./perm)
154          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155       elif perm.getRank()==2:              perm=perm*util.kronecker(self.domain.getDim())
156          perm_inv=util.inverse(perm)          
157              if self.useVPIteration:          
158                self.scale=1.           elif perm.getRank()==2:
159              else:              perm_inv=util.inverse(perm)
               self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())  
           perm_inv*=(1./self.scale)  
           perm=perm*self.scale  
      else:  
         raise ValueError,"illegal rank of permeability."  
   
      self.__permeability=perm  
      self.__permeability_inv=perm_inv  
   
      self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale  
          if self.useVPIteration:  
         if  self.solveForFlux:  
            self.__pde_k.setValue(D=self.__permeability_inv)  
         else:  
            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
         self.__pde_p.setValue(A=self.__permeability)  
160           else:           else:
161              D=self.__pde_k.createCoefficient("D")              raise ValueError("illegal rank of permeability.")
162              A=self.__pde_k.createCoefficient("A")          
163              D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv           self.__permeability=perm
164              for i in range(self.domain.getDim()):           self.__permeability_inv=perm_inv
165                 for j in range(self.domain.getDim()):      
166                   A[i,i,j,j]=self.__l2           #====================
167              A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability           self.__pde_p.setValue(A=self.__permeability)
168              self.__pde_k.setValue(A=A, D=D)           if self.solver  == self.EVAL:
169        if g !=None:                pass # no extra work required
170       g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))           elif self.solver  == self.POST:
171       if g.isEmpty():                k=util.kronecker(self.domain.getDim())
172            g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
173       else:                self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
174          if not g.getShape()==(self.domain.getDim(),):           elif self.solver  == self.SMOOTH:
175                raise ValueError,"illegal shape of g"              self.__pde_v.setValue(D=self.__permeability_inv)
176          self.__g=g  
177        elif permeability!=None:        if g != None:
178               X          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
179            if g.isEmpty():
180                 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
181            else:
182                 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
183            self.__g=g
184        if f !=None:        if f !=None:
185       f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
186       if f.isEmpty():           if f.isEmpty():      
187            f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
188       else:           else:
189           if f.getRank()>0: raise ValueError,"illegal rank of f."               if f.getRank()>0: raise ValueError("illegal rank of f.")
190           self.__f=f           self.__f=f
191    
192       def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          if self.__pde_v == None:
198              return None
199          else:
200              return self.__pde_v.getSolverOptions()
201                
202     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def setSolverOptionsFlux(self, options=None):
203        """        """
204        solves the problem.        Sets the solver options used to solve the flux problems
205          If ``options`` is not present, the options are reset to default
206          :param options: `SolverOptions`
207          """
208          if not self.__pde_v == None:
209              self.__pde_v.setSolverOptions(options)
210        
211       def getSolverOptionsPressure(self):
212          """
213          Returns the solver options used to solve the pressure problems
214          :return: `SolverOptions`
215          """
216          return self.__pde_p.getSolverOptions()
217                
218        The iteration is terminated if the residual norm is less then self.getTolerance().     def setSolverOptionsPressure(self, options=None):
   
       :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. `escript.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. `escript.Data`).  
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
       :return: flux and pressure  
       :rtype: ``tuple`` of `escript.Data`.  
   
       :note: The problem is solved as a least squares form  
              *(K^[-1]+D^* l2 D)u+G p=D^* l2 * 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}*.  
219        """        """
220        self.verbose=verbose        Sets the solver options used to solve the pressure problems
221        if self.useVPIteration:        If ``options`` is not present, the options are reset to default
222            return self.__solveVP(u0,p0,max_iter,max_num_corrections)        
223        else:        :param options: `SolverOptions`
224            X=self.__pde_k.createCoefficient("X")        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
225            Y=self.__pde_k.createCoefficient("Y")        """
226            Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)        return self.__pde_p.setSolverOptions(options)
227            rtmp=self.__f * self.__l2 * self.scale        
228            for i in range(self.domain.getDim()): X[i,i]=rtmp     def solve(self, u0, p0):
           X[self.domain.getDim(),:]=self.__g*self.scale  
           r=self.__pde_k.createCoefficient("r")  
           r[:self.domain.getDim()]=u0*self.scale  
           r[self.domain.getDim()]=p0  
           self.__pde_k.setValue(X=X, Y=Y, r=r)  
           self.__pde_k.getSolverOptions().setVerbosity(self.verbose)  
           #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)  
           self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)  
           U=self.__pde_k.getSolution()  
           # self.__pde_k.getOperator().saveMM("k.mm")  
           u=U[:self.domain.getDim()]*(1./self.scale)  
           p=U[self.domain.getDim()]  
           if self.verbose:  
         KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)  
         def_p=self.__g-(u+KGp)  
         def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))  
           return u,p  
   
    def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):  
229        """        """
230        solves the problem.        solves the problem.
231                
       The iteration is terminated if the residual norm is less than self.getTolerance().  
   
232        :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.
233        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
234        :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.
# Line 299  class DarcyFlow(object): Line 236  class DarcyFlow(object):
236        :return: flux and pressure        :return: flux and pressure
237        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
238    
       :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}*.  
239        """        """
240        rtol=self.getTolerance()        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
241        atol=self.getAbsoluteTolerance()                              Y=self.__f * 1./self.perm_scale,
242        self.setSubProblemTolerance()                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
243        num_corrections=0                              r=p0)
       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 __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l2  
       else:  
      return (f-util.div(v))*self.__l2  
    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)        
244        p=self.__pde_p.getSolution()        p=self.__pde_p.getSolution()
245        return p        u = self.getFlux(p, u0)
246          return u,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=escript.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=escript.Data())  
       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])  
         
    def __Aprod_p(self,dp):  
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())  
       du=self.__pde_k.getSolution()  
       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()    
247                
248        #v=self.__getFlux(p, u0, f=self.__f, g=g2)           def getFlux(self,p, u0=None):
    def __Msolve_PCG_p(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
       self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())  
       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,escript.Function(self.domain)))**2))  
   
    def __L2_r(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))  
   
    def setTolerance(self,rtol=1e-4):  
       """  
       sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
   
       *|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.):  
       """  
       sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if  
         
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
   
   
       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.)  
      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``    
         """  
         self.domain=domain  
         if weight == None:  
            s=self.domain.getSize()  
            self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2  
            # self.__l2=(3.*util.longestEdge(self.domain))**2  
            #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2  
         else:  
            self.__l2=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.__l2*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=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         self.__g=escript.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. `escript.Data`)  
         :param g: flux sources/sinks  
         :type g: vector values on the domain (e.g. `escript.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. `escript.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. `escript.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. `escript.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=escript.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=escript.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``  
249          """          """
250          if atol<0:          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
             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. `escript.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. `escript.Data`).  
          :param verbose: if set some information on iteration progress are printed  
          :type verbose: ``bool``  
          :return: flux and pressure  
          :rtype: ``tuple`` of `escript.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,escript.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,escript.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,escript.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=escript.Data(), r=escript.Data())  
           du=self.__pde_v.getSolution()  
           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=escript.Data(), r=escript.Data())  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=escript.Data()):  
         """  
         returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``  
251          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
252          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
253    
254          :param p: pressure.          :param p: pressure.
255          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
256          :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``.
257          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
258          :return: flux          :return: flux
259          :rtype: `escript.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}*  
260          """          """
261      self.setSubProblemTolerance()          if self.solver  == self.EVAL:
262          g=self.__g             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
263          f=self.__f          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
264          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
265          if p == None:              if u0 == None:
266             self.__pde_v.setValue(Y=g)                 self.__pde_v.setValue(r=escript.Data())
267          else:              else:
268             self.__pde_v.setValue(Y=g-self.__Q(p))                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
269          return self.__pde_v.getSolution()                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
270                   u= self.__pde_v.getSolution() * self.perm_scale
271            return u
272          
273  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
274       """       """
275       solves       solves
# Line 835  class StokesProblemCartesian(Homogeneous Line 288  class StokesProblemCartesian(Homogeneous
288              sp.setTolerance()              sp.setTolerance()
289              sp.initialize(...)              sp.initialize(...)
290              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
291                sp.setStokesEquation(...) # new values for some parameters
292                v1,p1=sp.solve(v,p)
293       """       """
294       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
295           """           """
# Line 849  class StokesProblemCartesian(Homogeneous Line 304  class StokesProblemCartesian(Homogeneous
304           """           """
305           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
306           self.domain=domain           self.domain=domain
307           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())
308           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
309            
310           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
311           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 858  class StokesProblemCartesian(Homogeneous Line 313  class StokesProblemCartesian(Homogeneous
313    
314           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
315           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
316       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
317           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
318    
319       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 867  class StokesProblemCartesian(Homogeneous Line 322  class StokesProblemCartesian(Homogeneous
322            
323       :rtype: `SolverOptions`       :rtype: `SolverOptions`
324       """       """
325       return self.__pde_u.getSolverOptions()           return self.__pde_v.getSolverOptions()
326       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
327           """           """
328       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 875  class StokesProblemCartesian(Homogeneous Line 330  class StokesProblemCartesian(Homogeneous
330       :param options: new solver  options       :param options: new solver  options
331       :type options: `SolverOptions`       :type options: `SolverOptions`
332       """       """
333           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
334       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
335           """           """
336       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
337       :rtype: `SolverOptions`       :rtype: `SolverOptions`
338       """       """
339       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
340       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
341           """           """
342       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
343       :param options: new solver  options       :param options: new solver  options
344       :type options: `SolverOptions`       :type options: `SolverOptions`
345       """       """
346       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
347    
348       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
349           """           """
# Line 898  class StokesProblemCartesian(Homogeneous Line 353  class StokesProblemCartesian(Homogeneous
353       :param options: new solver options       :param options: new solver options
354       :type options: `SolverOptions`       :type options: `SolverOptions`
355       """       """
356       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
357       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
358           """           """
359       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 906  class StokesProblemCartesian(Homogeneous Line 361  class StokesProblemCartesian(Homogeneous
361            
362       :rtype: `SolverOptions`       :rtype: `SolverOptions`
363       """       """
364       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
365    
366       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
367           """           """
368           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
369           :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.
370           """           """
371           pass           pass
372       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 933  class StokesProblemCartesian(Homogeneous Line 388  class StokesProblemCartesian(Homogeneous
388              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
389              kk=util.outer(k,k)              kk=util.outer(k,k)
390              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
391          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
392              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)))
393          if restoration_factor!=None:          if restoration_factor!=None:
394              n=self.domain.getNormal()              n=self.domain.getNormal()
395              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
396          if fixed_u_mask!=None:          if fixed_u_mask!=None:
397              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
398          if f!=None: self.__f=f          if f!=None: self.__f=f
399          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
400          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
# Line 957  class StokesProblemCartesian(Homogeneous Line 412  class StokesProblemCartesian(Homogeneous
412          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
413          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
414          :param stress: initial stress          :param stress: initial stress
415      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
416          """          """
417          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
418    
# Line 970  class StokesProblemCartesian(Homogeneous Line 425  class StokesProblemCartesian(Homogeneous
425           :rtype: ``float``           :rtype: ``float``
426           """           """
427           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
428       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
429       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
430           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
431           return out           return out
432    
# Line 1012  class StokesProblemCartesian(Homogeneous Line 467  class StokesProblemCartesian(Homogeneous
467    
468       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
469           """           """
470           return the value for v for a given p (overwrite)           return the value for v for a given p
471    
472           :param p: a pressure           :param p: a pressure
473           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
474           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
475           """           """
476           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
477           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
478       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
479       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
480           if self.__stress.isEmpty():           if self.__stress.isEmpty():
481              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)))
482           else:           else:
483              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)))
484           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
485           return  out           return  out
486    
487       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 1046  class StokesProblemCartesian(Homogeneous Line 501  class StokesProblemCartesian(Homogeneous
501           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
502           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
503           """           """
504           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
505           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
506           return  out           return  out
507    
508       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
# Line 1060  class StokesProblemCartesian(Homogeneous Line 515  class StokesProblemCartesian(Homogeneous
515           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
516           """           """
517           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
518       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
519       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
520           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
521           return out           return out

Legend:
Removed from v.3452  
changed lines
  Added in v.3852

  ViewVC Help
Powered by ViewVC 1.1.26