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

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

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

revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC revision 3462 by caltinay, Mon Feb 7 02:01:50 2011 UTC
# Line 32  Some models for flow Line 32  Some models for flow
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
# Line 49  class DarcyFlow(object): Line 49  class DarcyFlow(object):
49     :note: The problem is solved in a least squares formulation.     :note: The problem is solved in a least squares formulation.
50     """     """
51        
52     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):
53        """        """
54        initializes the Darcy flux problem        initializes the Darcy flux problem
55        :param domain: domain of the problem        :param domain: domain of the problem
# Line 60  class DarcyFlow(object): Line 60  class DarcyFlow(object):
60        :type adaptSubTolerance: ``bool``        :type adaptSubTolerance: ``bool``
61        :param solveForFlux: if True the solver solves for the flux (do not use!)        :param solveForFlux: if True the solver solves for the flux (do not use!)
62        :type solveForFlux: ``bool``          :type solveForFlux: ``bool``  
63          :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.
64          :type useVPIteration: ``bool``    
65        """        """
66        self.domain=domain        self.domain=domain
67        self.solveForFlux=solveForFlux        self.useVPIteration=useVPIteration
68        self.useReduced=useReduced        self.useReduced=useReduced
69        self.__adaptSubTolerance=adaptSubTolerance        self.weighting_scale=weighting_scale
70        self.verbose=False        if self.useVPIteration:
71                   self.solveForFlux=solveForFlux
72        self.__pde_k=LinearPDESystem(domain)           self.__adaptSubTolerance=adaptSubTolerance
73        self.__pde_k.setSymmetryOn()           self.verbose=False
74        if self.useReduced: self.__pde_k.setReducedOrderOn()          
75             self.__pde_k=LinearPDESystem(domain)
76        self.__pde_p=LinearSinglePDE(domain)           self.__pde_k.setSymmetryOn()
77        self.__pde_p.setSymmetryOn()           if self.useReduced: self.__pde_k.setReducedOrderOn()
78        if self.useReduced: self.__pde_p.setReducedOrderOn()    
79             self.__pde_p=LinearSinglePDE(domain)
80        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting           self.__pde_p.setSymmetryOn()
81        # self.__pde_l.setSymmetryOn()           if self.useReduced: self.__pde_p.setReducedOrderOn()
82        # if self.useReduced: self.__pde_l.setReducedOrderOn()           self.setTolerance()
83        self.setTolerance()           self.setAbsoluteTolerance()
84        self.setAbsoluteTolerance()        else:
85        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))           self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
86        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))           self.__pde_k.setSymmetryOn()
87             if self.useReduced: self.__pde_k.setReducedOrderOn()
88             C=self.__pde_k.createCoefficient("C")
89             B=self.__pde_k.createCoefficient("B")
90             for i in range(self.domain.getDim()):
91                C[i,self.domain.getDim(),i]=1
92                B[self.domain.getDim(),i,i]=1
93             self.__pde_k.setValue(C=C, B=B)
94          self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
95          self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
96                
97     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
98        """        """
# Line 129  class DarcyFlow(object): Line 140  class DarcyFlow(object):
140        """        """
141        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
142    
    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)  
   
143    
144     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):
145        """        """
146        assigns values to model parameters        assigns values to model parameters
147          
148        :param f: volumetic sources/sinks        :param f: volumetic sources/sinks
149        :type f: scalar value on the domain (e.g. `Data`)        :type f: scalar value on the domain (e.g. `escript.Data`)
150        :param g: flux sources/sinks        :param g: flux sources/sinks
151        :type g: vector values on the domain (e.g. `Data`)        :type g: vector values on the domain (e.g. `escript.Data`)
152        :param location_of_fixed_pressure: mask for locations where pressure is fixed        :param location_of_fixed_pressure: mask for locations where pressure is fixed
153        :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)        :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
154        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
155        :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
156        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
157        ``s`` on the main diagonal is used.        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
158        :type permeability: scalar or tensor values on the domain (e.g. `Data`)  
159        :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.
160        :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)        :note: at any point on the boundary of the domain the pressure
161        or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal               (``location_of_fixed_pressure`` >0) or the normal component of the
162        is along the *x_i* axis.               flux (``location_of_fixed_flux[i]>0``) if direction of the normal
163                 is along the *x_i* axis.
164    
165        """        """
166        if f !=None:        if self.useVPIteration:
167       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
168       if f.isEmpty():           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)
169          f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        else:
170       else:           q=self.__pde_k.getCoefficient("q")
171          if f.getRank()>0: raise ValueError,"illegal rank of f."           if q.isEmpty(): q=self.__pde_k.createCoefficient("q")
172          self.__f=f           if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure
173        if g !=None:           if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux
174       g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           self.__pde_k.setValue(q=q)
      if g.isEmpty():  
         g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
            raise ValueError,"illegal shape of g"  
         self.__g=g  
       if location_of_fixed_pressure!=None:  
            self.__pde_p.setValue(q=location_of_fixed_pressure)  
            #self.__pde_l.setValue(q=location_of_fixed_pressure)  
       if location_of_fixed_flux!=None:  
            self.__pde_k.setValue(q=location_of_fixed_flux)  
175                            
176          # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)
177        if permeability!=None:        if permeability!=None:
178       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))
179             V=util.vol(self.domain)
180       if perm.getRank()==0:       if perm.getRank()==0:
181          perm_inv=(1./perm)*util.kronecker(self.domain.getDim())          perm_inv=(1./perm)
182          perm=perm*util.kronecker(self.domain.getDim())              if self.useVPIteration:
183                  self.scale=1.
184                else:
185                  self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)
186    
187            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
188            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
189       elif perm.getRank()==2:       elif perm.getRank()==2:
190          perm_inv=util.inverse(perm)          perm_inv=util.inverse(perm)
191                if self.useVPIteration:
192                  self.scale=1.
193                else:
194                  self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())
195              perm_inv*=(1./self.scale)
196              perm=perm*self.scale
197       else:       else:
198          raise ValueError,"illegal rank of permeability."          raise ValueError,"illegal rank of permeability."
199    
200       self.__permeability=perm       self.__permeability=perm
201       self.__permeability_inv=perm_inv       self.__permeability_inv=perm_inv
      self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10  
      #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10  
      if  self.solveForFlux:  
         self.__pde_k.setValue(D=self.__permeability_inv)  
      else:  
         self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
      self.__pde_p.setValue(A=self.__permeability)  
      #self.__pde_l.setValue(D=1/self.__l)  
          #self.__pde_l.setValue(A=self.__permeability)  
   
    def __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l  
       else:  
      return (f-util.div(v))*self.__l  
       # if f == None:  
       #      self.__pde_l.setValue(Y=-util.div(v))    
       # else:  
       #      return (f-util.div(v))/self.__l  
       # return self.__pde_l.getSolution()  
         
    def __getPressure(self, v, p0, g=None):  
       # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0  
       if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"  
       if g == None:  
      self.__pde_p.setValue(X=-v, r=p0)  
       else:  
      self.__pde_p.setValue(X=g-v, r=p0)        
       p=self.__pde_p.getSolution()  
       return p  
   
    def __Aprod_v(self,dv):  
       # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv    
       dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)  
       a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp  
       b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)  
       return ArithmeticTuple(a,b,-dp)  
202    
203     def __Msolve_PCG_v(self,r):       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
204        # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv           if self.useVPIteration:
205        if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"          if  self.solveForFlux:
206        self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())             self.__pde_k.setValue(D=self.__permeability_inv)
207        # self.__pde_p.getOperator().saveMM("prec.mm")          else:
208        return self.__pde_k.getSolution()             self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
209            self.__pde_p.setValue(A=self.__permeability)
210             else:
211                D=self.__pde_k.createCoefficient("D")
212                A=self.__pde_k.createCoefficient("A")
213                D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
214                for i in range(self.domain.getDim()):
215                   for j in range(self.domain.getDim()):
216                     A[i,i,j,j]=self.__l2
217                A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
218                self.__pde_k.setValue(A=A, D=D)
219          if g !=None:
220         g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))
221         if g.isEmpty():
222              g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))
223         else:
224            if not g.getShape()==(self.domain.getDim(),):
225                  raise ValueError,"illegal shape of g"
226            self.__g=g
227          elif permeability!=None:
228                 X
229          if f !=None:
230         f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))
231         if f.isEmpty():
232              f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))
233         else:
234             if f.getRank()>0: raise ValueError,"illegal rank of f."
235             self.__f=f
236    
    def __inner_PCG_v(self,v,r):  
       return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])  
237                
    def __Aprod_p(self,dp):  
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())  
       du=self.__pde_k.getSolution()  
       # self.__pde_v.getOperator().saveMM("proj.mm")  
       return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)  
   
    def __getFlux(self,p, v0, f=None, g=None):  
       # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp  
       if f!=None:  
      self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))  
       self.__pde_k.setValue(r=v0)  
       g2=util.tensor_mult(self.__permeability_inv,g)  
       if p == None:  
      self.__pde_k.setValue(Y=g2)  
       else:  
      self.__pde_k.setValue(Y=g2-util.grad(p))  
       return self.__pde_k.getSolution()    
         
       #v=self.__getFlux(p, u0, f=self.__f, g=g2)        
    def __Msolve_PCG_p(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
       self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
   
    def __L2(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
238     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
239        """        """
240        solves the problem.        solves the problem.
# Line 296  class DarcyFlow(object): Line 242  class DarcyFlow(object):
242        The iteration is terminated if the residual norm is less then self.getTolerance().        The iteration is terminated if the residual norm is less then self.getTolerance().
243    
244        :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.
245        :type u0: vector value on the domain (e.g. `Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
246        :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.
247        :type p0: scalar value on the domain (e.g. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
248        :param verbose: if set some information on iteration progress are printed        :param verbose: if set some information on iteration progress are printed
249        :type verbose: ``bool``        :type verbose: ``bool``
250        :return: flux and pressure        :return: flux and pressure
251        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.Data`.
252    
253        :note: The problem is solved as a least squares form        :note: The problem is solved as a least squares form
254        *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*               *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
255        *G^*u+G^* K Gp=G^*g*               *G^*u+*G^* K Gp=G^*g*
256                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
       where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.  
257        """        """
258        self.verbose=verbose        self.verbose=verbose
259          if self.useVPIteration:
260              return self.__solveVP(u0,p0,max_iter,max_num_corrections)
261          else:
262              X=self.__pde_k.createCoefficient("X")
263              Y=self.__pde_k.createCoefficient("Y")
264              Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)
265              rtmp=self.__f * self.__l2 * self.scale
266              for i in range(self.domain.getDim()): X[i,i]=rtmp
267              X[self.domain.getDim(),:]=self.__g*self.scale
268              r=self.__pde_k.createCoefficient("r")
269              r[:self.domain.getDim()]=u0*self.scale
270              r[self.domain.getDim()]=p0
271              self.__pde_k.setValue(X=X, Y=Y, r=r)
272              self.__pde_k.getSolverOptions().setVerbosity(self.verbose)
273              #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)
274              self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)
275              U=self.__pde_k.getSolution()
276              # self.__pde_k.getOperator().saveMM("k.mm")
277              u=U[:self.domain.getDim()]*(1./self.scale)
278              p=U[self.domain.getDim()]
279              if self.verbose:
280            KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)
281            def_p=self.__g-(u+KGp)
282            def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))
283            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))
284            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
285              return u,p
286    
287       def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
288          """
289          solves the problem.
290          
291          The iteration is terminated if the residual norm is less than self.getTolerance().
292    
293          :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.
294          :type u0: vector value on the domain (e.g. `escript.Data`).
295          :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.
296          :type p0: scalar value on the domain (e.g. `escript.Data`).
297          :return: flux and pressure
298          :rtype: ``tuple`` of `escript.Data`.
299    
300          :note: The problem is solved as a least squares form
301                 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
302                 *G^*u+*G^* K Gp=G^*g*
303                 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
304          """
305        rtol=self.getTolerance()        rtol=self.getTolerance()
306        atol=self.getAbsoluteTolerance()        atol=self.getAbsoluteTolerance()
307        self.setSubProblemTolerance()        self.setSubProblemTolerance()
# Line 372  class DarcyFlow(object): Line 363  class DarcyFlow(object):
363          if self.verbose: print "DarcyFlux: stopping criterium reached."          if self.verbose: print "DarcyFlux: stopping criterium reached."
364          converged=True          converged=True
365        return v,p+p0        return v,p+p0
366    
367       def __applWeight(self, v, f=None):
368          # solves L p = f-Dv with p = 0
369          if self.verbose: print "DarcyFlux: Applying weighting operator"
370          if f == None:
371         return -util.div(v)*self.__l2
372          else:
373         return (f-util.div(v))*self.__l2
374       def __getPressure(self, v, p0, g=None):
375          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
376          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
377          if g == None:
378         self.__pde_p.setValue(X=-v, r=p0)
379          else:
380         self.__pde_p.setValue(X=g-v, r=p0)      
381          p=self.__pde_p.getSolution()
382          return p
383    
384       def __Aprod_v(self,dv):
385          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
386          dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)
387          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
388          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
389          return ArithmeticTuple(a,b,-dp)
390    
391       def __Msolve_PCG_v(self,r):
392          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
393          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
394          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())
395          return self.__pde_k.getSolution()
396    
397       def __inner_PCG_v(self,v,r):
398          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
399          
400       def __Aprod_p(self,dp):
401          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
402          Gdp=util.grad(dp)
403          self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())
404          du=self.__pde_k.getSolution()
405          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
406    
407       def __getFlux(self,p, v0, f=None, g=None):
408          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
409          if f!=None:
410         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
411          self.__pde_k.setValue(r=v0)
412          g2=util.tensor_mult(self.__permeability_inv,g)
413          if p == None:
414         self.__pde_k.setValue(Y=g2)
415          else:
416         self.__pde_k.setValue(Y=g2-util.grad(p))
417          return self.__pde_k.getSolution()  
418          
419          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
420       def __Msolve_PCG_p(self,r):
421          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
422          self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data())
423          return self.__pde_p.getSolution()
424            
425       def __inner_PCG_p(self,p,r):
426           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
427    
428       def __L2(self,v):
429          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
430    
431       def __L2_r(self,v):
432          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))
433    
434     def setTolerance(self,rtol=1e-4):     def setTolerance(self,rtol=1e-4):
435        """        """
436        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
# Line 433  class DarcyFlow(object): Line 492  class DarcyFlow(object):
492       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
493       self.getSolverOptionsPressure().setTolerance(sub_tol)       self.getSolverOptionsPressure().setTolerance(sub_tol)
494       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
      self.getSolverOptionsWeighting().setTolerance(sub_tol)  
      self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)  
495       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
496    
497    
# Line 463  class DarcyFlowOld(object): Line 520  class DarcyFlowOld(object):
520          self.domain=domain          self.domain=domain
521          if weight == None:          if weight == None:
522             s=self.domain.getSize()             s=self.domain.getSize()
523             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
524             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l2=(3.*util.longestEdge(self.domain))**2
525             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
526          else:          else:
527             self.__l=weight             self.__l2=weight
528          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
529          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
530          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
531          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain)))
532          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
533          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
534          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
535          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
536          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
537          self.setTolerance()          self.setTolerance()
538          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
539      self.__adaptSubTolerance=adaptSubTolerance      self.__adaptSubTolerance=adaptSubTolerance
# Line 527  class DarcyFlowOld(object): Line 584  class DarcyFlowOld(object):
584          assigns values to model parameters          assigns values to model parameters
585    
586          :param f: volumetic sources/sinks          :param f: volumetic sources/sinks
587          :type f: scalar value on the domain (e.g. `Data`)          :type f: scalar value on the domain (e.g. `escript.Data`)
588          :param g: flux sources/sinks          :param g: flux sources/sinks
589          :type g: vector values on the domain (e.g. `Data`)          :type g: vector values on the domain (e.g. `escript.Data`)
590          :param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
591          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
592          :param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
593          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
594          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
595                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
596                               ``v`` on the main diagonal is used.                               ``v`` on the main diagonal is used.
597          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)          :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`)
598    
599          :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.
600          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
# Line 547  class DarcyFlowOld(object): Line 604  class DarcyFlowOld(object):
604          if f !=None:          if f !=None:
605             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
606             if f.isEmpty():             if f.isEmpty():
607                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
608             else:             else:
609                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
610             self.__f=f             self.__f=f
611          if g !=None:          if g !=None:
612             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
613             if g.isEmpty():             if g.isEmpty():
614               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
615             else:             else:
616               if not g.getShape()==(self.domain.getDim(),):               if not g.getShape()==(self.domain.getDim(),):
617                 raise ValueError,"illegal shape of g"                 raise ValueError,"illegal shape of g"
# Line 647  class DarcyFlowOld(object): Line 704  class DarcyFlowOld(object):
704           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
705    
706           :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.
707           :type u0: vector value on the domain (e.g. `Data`).           :type u0: vector value on the domain (e.g. `escript.Data`).
708           :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.
709           :type p0: scalar value on the domain (e.g. `Data`).           :type p0: scalar value on the domain (e.g. `escript.Data`).
710           :param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
711           :type verbose: ``bool``           :type verbose: ``bool``
712           :return: flux and pressure           :return: flux and pressure
713           :rtype: ``tuple`` of `Data`.           :rtype: ``tuple`` of `escript.Data`.
714    
715           :note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
716    
# Line 699  class DarcyFlowOld(object): Line 756  class DarcyFlowOld(object):
756                 if self.verbose:                 if self.verbose:
757                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
758                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
759                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),)
760                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
761                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
762                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
763                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
764                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
765                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     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)
766                     num_corrections+=1                     num_corrections+=1
767                 else:                 else:
768                     converged=True                     converged=True
769           return v,p           return v,p
770      def __L2(self,v):      def __L2(self,v):
771           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
772    
773      def __Q(self,p):      def __Q(self,p):
774            return util.tensor_mult(self.__permeability,util.grad(p))            return util.tensor_mult(self.__permeability,util.grad(p))
# Line 719  class DarcyFlowOld(object): Line 776  class DarcyFlowOld(object):
776      def __Aprod(self,dp):      def __Aprod(self,dp):
777            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
778            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
779            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data())
780            du=self.__pde_v.getSolution()            du=self.__pde_v.getSolution()
           # self.__pde_v.getOperator().saveMM("proj.mm")  
781            return Qdp+du            return Qdp+du
782      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
783           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 731  class DarcyFlowOld(object): Line 787  class DarcyFlowOld(object):
787    
788      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
789        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
790            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            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")  
791            return self.__pde_p.getSolution()            return self.__pde_p.getSolution()
792    
793      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
794          """          """
795          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 ``fixed_flux``
796          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
797          Note that ``g`` and ``f`` are used, see `setValue`.          Note that ``g`` and ``f`` are used, see `setValue`.
798    
799          :param p: pressure.          :param p: pressure.
800          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
801          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
802          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).
803          :return: flux          :return: flux
804          :rtype: `Data`          :rtype: `escript.Data`
805          :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}*          :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}*
806                 for the permeability *k_{ij}*                 for the permeability *k_{ij}*
807          """          """
808      self.setSubProblemTolerance()      self.setSubProblemTolerance()
809          g=self.__g          g=self.__g
810          f=self.__f          f=self.__f
811          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)
812          if p == None:          if p == None:
813             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
814          else:          else:
# Line 875  class StokesProblemCartesian(Homogeneous Line 930  class StokesProblemCartesian(Homogeneous
930          if eta !=None:          if eta !=None:
931              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
932              kk=util.outer(k,k)              kk=util.outer(k,k)
933              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
934          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
935              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
936          if restoration_factor!=None:          if restoration_factor!=None:
# Line 887  class StokesProblemCartesian(Homogeneous Line 942  class StokesProblemCartesian(Homogeneous
942          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
943          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
944    
945       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
946          """          """
947          assigns values to the model parameters          assigns values to the model parameters
948    
# Line 927  class StokesProblemCartesian(Homogeneous Line 982  class StokesProblemCartesian(Homogeneous
982           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
983           :rtype: ``float``           :rtype: ``float``
984           """           """
985           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
986    
987       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
988           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 993  class StokesProblemCartesian(Homogeneous
993           :return: inner product of p0 and p1           :return: inner product of p0 and p1
994           :rtype: ``float``           :rtype: ``float``
995           """           """
996           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
997           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
998           return util.integrate(s0*s1)           return util.integrate(s0*s1)
999    
1000       def norm_v(self,v):       def norm_v(self,v):
# Line 979  class StokesProblemCartesian(Homogeneous Line 1034  class StokesProblemCartesian(Homogeneous
1034          :rtype: equal to the type of p          :rtype: equal to the type of p
1035          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1036          """          """
1037          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
1038    
1039       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
1040           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 1044  class StokesProblemCartesian(Homogeneous
1044           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
1045           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1046           """           """
1047           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
1048           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
1049           return  out           return  out
1050    

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

  ViewVC Help
Powered by ViewVC 1.1.26