/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3432 by jfenwick, Fri Jan 7 01:32:07 2011 UTC revision 3515 by gross, Thu May 19 08:20:57 2011 UTC
# 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 SIMPLE: simple solver
50       :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52       :cvar SYMSTAB: solver uses symmetric stabilization
53     """     """
54         SIMPLE="SIMPLE"
55     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     POST="POST"
56       STAB="STAB"
57       SYMSTAB="SYMSTAB"
58       def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59        """        """
60        initializes the Darcy flux problem        initializes the Darcy flux problem
61        :param domain: domain of the problem        :param domain: domain of the problem
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.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
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        self.domain=domain        self.domain=domain
74        self.solveForFlux=solveForFlux        self.solver=solver
75        self.useReduced=useReduced        self.useReduced=useReduced
76        self.__adaptSubTolerance=adaptSubTolerance        self.verbose=verbose
77        self.verbose=False        self.scale=1.
78                
79        self.__pde_k=LinearPDESystem(domain)        
80        self.__pde_k.setSymmetryOn()        self.__pde_v=LinearPDESystem(domain)
81        if self.useReduced: self.__pde_k.setReducedOrderOn()        self.__pde_v.setSymmetryOn()
82          if self.useReduced: self.__pde_v.setReducedOrderOn()
83        self.__pde_p=LinearSinglePDE(domain)        self.__pde_p=LinearSinglePDE(domain)
84        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
85        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
86          
87        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting        if self.solver  == self.SIMPLE:
88        # self.__pde_l.setSymmetryOn()       if self.verbose: print "DarcyFlow: simple solver is used."
89        # if self.useReduced: self.__pde_l.setReducedOrderOn()           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90        self.setTolerance()        elif self.solver  == self.POST:
91        self.setAbsoluteTolerance()       self.w=w
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          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104          self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105          self.setTolerance()
106        
107            
108       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109          """
110          assigns values to model parameters
111    
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122    
123          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136                
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140        
141         perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142             V=util.vol(self.domain)
143             l=V**(1./self.domain.getDim())
144            
145         if perm.getRank()==0:
146            perm_inv=(1./perm)
147                self.scale=util.integrate(perm_inv)/V*l
148            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150            
151            
152         elif perm.getRank()==2:
153            perm_inv=util.inverse(perm)
154                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155            perm_inv*=(1./self.scale)
156            perm=perm*self.scale
157         else:
158            raise ValueError,"illegal rank of permeability."
159            
160         self.__permeability=perm
161         self.__permeability_inv=perm_inv
162         if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163        
164         if self.solver  == self.SIMPLE:
165            self.__pde_p.setValue(A=self.__permeability)
166         elif self.solver  == self.POST:
167            self.__pde_p.setValue(A=self.__permeability)
168            k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172            self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():      
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
193        """        """
194        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
195        :return: `SolverOptions`        :return: `SolverOptions`
196        """        """
197        return self.__pde_k.getSolverOptions()        return self.__pde_v.getSolverOptions()
198                
199     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
200        """        """
201        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
202        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
203        :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.  
204        """        """
205        return self.__pde_v.setSolverOptions(options)        return 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=escript.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=escript.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=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())  
       # 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])  
         
    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()  
       # 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()    
223                
       #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=escript.Data(), r=escript.Data(), y=escript.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,escript.Function(self.domain)))**2))  
   
    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  
       *(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  
224     def setTolerance(self,rtol=1e-4):     def setTolerance(self,rtol=1e-4):
225        """        """
226        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
   
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
         
       where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).  
227                
228        :param rtol: relative tolerance for the pressure        :param rtol: relative tolerance for the pressure
229        :type rtol: non-negative ``float``        :type rtol: non-negative ``float``
# Line 385  class DarcyFlow(object): Line 231  class DarcyFlow(object):
231        if rtol<0:        if rtol<0:
232       raise ValueError,"Relative tolerance needs to be non-negative."       raise ValueError,"Relative tolerance needs to be non-negative."
233        self.__rtol=rtol        self.__rtol=rtol
234          
235     def getTolerance(self):     def getTolerance(self):
236        """        """
237        returns the relative tolerance        returns the relative tolerance
# Line 392  class DarcyFlow(object): Line 239  class DarcyFlow(object):
239        :rtype: ``float``        :rtype: ``float``
240        """        """
241        return self.__rtol        return self.__rtol
242          
243     def setAbsoluteTolerance(self,atol=0.):     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244        """        """
245        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        solves the problem.
246                
247        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*        The iteration is terminated if the residual norm is less then self.getTolerance().
   
   
       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}*.  
248    
249        :param atol: absolute tolerance for the pressure        :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.
250        :type atol: non-negative ``float``        :type u0: vector value on the domain (e.g. `escript.Data`).
251        """        :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.
252        if atol<0:        :type p0: scalar value on the domain (e.g. `escript.Data`).
253       raise ValueError,"Absolute tolerance needs to be non-negative."        :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254        self.__atol=atol        :type max_iter: ``int``
255     def getAbsoluteTolerance(self):        :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256        """                             A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257        returns the absolute tolerance                             of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258        :return: current absolute tolerance        :type iter_restart: ``int``
259        :rtype: ``float``        :return: flux and pressure
260        """        :rtype: ``tuple`` of `escript.Data`.
       return self.__atol  
    def getSubProblemTolerance(self):  
       """  
       Returns a suitable subtolerance  
       :type: ``float``  
       """  
       return max(util.EPSILON**(0.5),self.getTolerance()**2)  
261    
    def setSubProblemTolerance(self):  
       """  
       Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
262        """        """
263        if self.__adaptSubTolerance:        # rescale initial guess:
264       sub_tol=self.getSubProblemTolerance()        p0=p0/self.scale
265       self.getSolverOptionsFlux().setTolerance(sub_tol)        if self.solver  == self.SIMPLE or self.solver  == self.POST :
266       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)          self.__pde_p.setValue(X=self.__g ,
267       self.getSolverOptionsPressure().setTolerance(sub_tol)                                Y=self.__f,
268       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269       self.getSolverOptionsWeighting().setTolerance(sub_tol)                                r=p0)
270       self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)          p=self.__pde_p.getSolution()
271       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol          u = self.getFlux(p, u0)
272          elif  self.solver  == self.STAB:
273        u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274  class DarcyFlowOld(object):        elif  self.solver  == self.SYMSTAB:
275      """      u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
     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.__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=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  
276            
277      *(I+D^*D)u=F*        if self.verbose:
278                KGp=util.tensor_mult(self.__permeability,util.grad(p))
279      If ``options`` is not present, the options are reset to default          def_p=self.__g-(u+KGp)
280      :param options: `SolverOptions`          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282      """          print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283      return self.__pde_v.setSolverOptions(options)        #rescale result
284      def getSolverOptionsPressure(self):        p=p*self.scale
285      """        return u,p
286      Returns the solver options used to solve the pressure problems        
287           def getFlux(self,p, u0=None):
     *(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=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``  
         """  
         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,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()  
           # 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=escript.Data(), r=escript.Data())  
           # self.__pde_p.getOperator().saveMM("prec.mm")  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=escript.Data()):  
288          """          """
289          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` and ``f`` are used, see `setValue`.
292    
293          :param p: pressure.          :param p: pressure.
294          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
295          :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``.
296          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          :return: flux          :return: flux
298          :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}*  
299          """          """
300      self.setSubProblemTolerance()          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301          g=self.__g              KGp=util.tensor_mult(self.__permeability,util.grad(p))
302          f=self.__f              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)              if u0 == None:
304          if p == None:             self.__pde_v.setValue(r=escript.Data())
305             self.__pde_v.setValue(Y=g)          else:
306          else:             self.__pde_v.setValue(r=u0)
307             self.__pde_v.setValue(Y=g-self.__Q(p))              u= self.__pde_v.getSolution()
308          return self.__pde_v.getSolution()      elif self.solver  == self.POST:
309                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311                if u0 == None:
312               self.__pde_v.setValue(r=escript.Data())
313            else:
314               self.__pde_v.setValue(r=u0)
315                u= self.__pde_v.getSolution()
316        elif self.solver  == self.STAB:
317             gp=util.grad(p)
318             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319                                   X= p * util.kronecker(self.domain.getDim()),
320                                   y= - p * self.domain.getNormal())                          
321             if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323             else:
324               self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326        elif  self.solver  == self.SYMSTAB:
327             gp=util.grad(p)
328             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330                                   y= escript.Data() )                          
331             if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333             else:
334               self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336        return u
337          
338        
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)  
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          p=GMRES(dp,
371                  self.__SYMSTAB_Aprod,
372              p0,
373              self.__inner,
374              atol=self.__norm(p0+dp)*self.getTolerance() ,
375              rtol=0.,
376              iter_max=max_iter,
377              iter_restart=iter_restart,
378              verbose=self.verbose,P_R=None)
379                
380              u=self.getFlux(p, u0)
381              return u,p
382    
383       def __L2(self,v):
384             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
385      
386       def __norm(self,r):
387             return util.sqrt(self.__inner(r,r))
388            
389       def __inner(self,r,s):
390             return util.integrate(util.inner(r,s), escript.Function(self.domain))
391            
392       def __STAB_Aprod(self,p):
393          gp=util.grad(p)
394          self.__pde_v.setValue(Y=-0.5*gp,
395                                X=-p*util.kronecker(self.__pde_v.getDomain()),
396                                y= p * self.domain.getNormal(),  
397                                r=escript.Data())
398          u = -self.__pde_v.getSolution()
399          self.__pde_p.setValue(Y=util.div(u),
400                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401                                y=escript.Data(),
402                                r=escript.Data())
403        
404          return  self.__pde_p.getSolution()
405      
406       def __SYMSTAB_Aprod(self,p):
407          gp=util.grad(p)
408          self.__pde_v.setValue(Y=0.5*gp ,
409                                X=escript.Data(),
410                                y=escript.Data(),  
411                                r=escript.Data())
412          u = -self.__pde_v.getSolution()
413          self.__pde_p.setValue(Y=escript.Data(),
414                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415                                y=escript.Data(),
416                                r=escript.Data())
417        
418          return  self.__pde_p.getSolution()
419          
420    
421  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
422       """       """
# Line 791  class StokesProblemCartesian(Homogeneous Line 450  class StokesProblemCartesian(Homogeneous
450           """           """
451           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
452           self.domain=domain           self.domain=domain
453           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())
454           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
455            
456           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
457           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 809  class StokesProblemCartesian(Homogeneous Line 468  class StokesProblemCartesian(Homogeneous
468            
469       :rtype: `SolverOptions`       :rtype: `SolverOptions`
470       """       """
471       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
472       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
473           """           """
474       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 817  class StokesProblemCartesian(Homogeneous Line 476  class StokesProblemCartesian(Homogeneous
476       :param options: new solver  options       :param options: new solver  options
477       :type options: `SolverOptions`       :type options: `SolverOptions`
478       """       """
479           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
480       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
481           """           """
482       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 876  class StokesProblemCartesian(Homogeneous Line 535  class StokesProblemCartesian(Homogeneous
535              kk=util.outer(k,k)              kk=util.outer(k,k)
536              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
537          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
538              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)))
539          if restoration_factor!=None:          if restoration_factor!=None:
540              n=self.domain.getNormal()              n=self.domain.getNormal()
541              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542          if fixed_u_mask!=None:          if fixed_u_mask!=None:
543              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
544          if f!=None: self.__f=f          if f!=None: self.__f=f
545          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
546          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
# Line 961  class StokesProblemCartesian(Homogeneous Line 620  class StokesProblemCartesian(Homogeneous
620           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
621           """           """
622           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
623           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
625       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626           if self.__stress.isEmpty():           if self.__stress.isEmpty():
627              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)))
628           else:           else:
629              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)))
630           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
631           return  out           return  out
632    
633       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 988  class StokesProblemCartesian(Homogeneous Line 647  class StokesProblemCartesian(Homogeneous
647           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
648           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
649           """           """
650           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))
651           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
652           return  out           return  out
653    
654       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

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

  ViewVC Help
Powered by ViewVC 1.1.26