/[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 3432 by jfenwick, Fri Jan 7 01:32:07 2011 UTC revision 3440 by gross, Fri Jan 14 00:04:53 2011 UTC
# Line 37  import util Line 37  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=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))           self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)
88        self.__g=escript.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 ``s`` on the main diagonal is used.        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
159        :type permeability: scalar or tensor values on the domain (e.g. `Data`)        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)
160    
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 (``location_of_fixed_pressure`` >0)
163        or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal is along the *x_i* axis.        or the normal component of the flux (``location_of_fixed_flux[i]>0``) if direction of the normal 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=escript.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=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
            raise ValueError,"illegal shape of g"  
         self.__g=g  
       if location_of_fixed_pressure!=None:  
            self.__pde_p.setValue(q=location_of_fixed_pressure)  
            #self.__pde_l.setValue(q=location_of_fixed_pressure)  
       if location_of_fixed_flux!=None:  
            self.__pde_k.setValue(q=location_of_fixed_flux)  
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
202       self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10  
203       #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
204       if  self.solveForFlux:           if self.useVPIteration:
205          self.__pde_k.setValue(D=self.__permeability_inv)          if  self.solveForFlux:
206               self.__pde_k.setValue(D=self.__permeability_inv)
207            else:
208               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:       else:
234          self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))           if f.getRank()>0: raise ValueError,"illegal rank of f."
235       self.__pde_p.setValue(A=self.__permeability)           self.__f=f
      #self.__pde_l.setValue(D=1/self.__l)  
          #self.__pde_l.setValue(A=self.__permeability)  
236    
    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()  
237                
238     def __getPressure(self, v, p0, g=None):     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
239        # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0        """
240        if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"        solves the problem.
       if g == None:  
      self.__pde_p.setValue(X=-v, r=p0)  
       else:  
      self.__pde_p.setValue(X=g-v, r=p0)        
       p=self.__pde_p.getSolution()  
       return p  
   
    def __Aprod_v(self,dv):  
       # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv    
       dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)  
       a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp  
       b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)  
       return ArithmeticTuple(a,b,-dp)  
   
    def __Msolve_PCG_v(self,r):  
       # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv  
       if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"  
       self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_k.getSolution()  
   
    def __inner_PCG_v(self,v,r):  
       return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])  
241                
242     def __Aprod_p(self,dp):        The iteration is terminated if the residual norm is less then self.getTolerance().
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())  
       du=self.__pde_k.getSolution()  
       # self.__pde_v.getOperator().saveMM("proj.mm")  
       return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)  
243    
244     def __getFlux(self,p, v0, f=None, g=None):        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
245        # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp        :type u0: vector value on the domain (e.g. `escript.Data`).
246        if f!=None:        :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       self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))        :type p0: scalar value on the domain (e.g. `escript.Data`).
248        self.__pde_k.setValue(r=v0)        :param verbose: if set some information on iteration progress are printed
249        g2=util.tensor_mult(self.__permeability_inv,g)        :type verbose: ``bool``
250        if p == None:        :return: flux and pressure
251       self.__pde_k.setValue(Y=g2)        :rtype: ``tuple`` of `escript.Data`.
       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=escript.Data(), r=escript.Data(), y=escript.Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
252    
253     def __L2(self,v):        :note: The problem is solved as a least squares form
254        return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))        *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g*
255          *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}*.
257          """
258          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 solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):
288        """        """
289        solves the problem.        solves the problem.
290                
291        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().
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.        :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. `Data`).        :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.        :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. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
297        :param verbose: if set some information on iteration progress are printed        :param verbose: if set some information on iteration progress are printed
298        :type verbose: ``bool``        :type verbose: ``bool``
299        :return: flux and pressure        :return: flux and pressure
300        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.Data`.
301    
302        :note: The problem is solved as a least squares form        :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*        *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
304        *G^*u+*G^* K Gp=G^*g*        *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}*.        where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
306        """        """
       self.verbose=verbose  
307        rtol=self.getTolerance()        rtol=self.getTolerance()
308        atol=self.getAbsoluteTolerance()        atol=self.getAbsoluteTolerance()
309        self.setSubProblemTolerance()        self.setSubProblemTolerance()
# Line 371  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 432  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 462  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()
# Line 526  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 646  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 720  class DarcyFlowOld(object): Line 780  class DarcyFlowOld(object):
780            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
781            self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.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 790  class DarcyFlowOld(object):
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=escript.Data(), r=escript.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=escript.Data()):      def getFlux(self,p=None, fixed_flux=escript.Data()):
# Line 741  class DarcyFlowOld(object): Line 799  class DarcyFlowOld(object):
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:

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

  ViewVC Help
Powered by ViewVC 1.1.26