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

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

  ViewVC Help
Powered by ViewVC 1.1.26