/[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 3452 by caltinay, Tue Jan 25 01:53:57 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
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.               (``location_of_fixed_pressure`` >0) or the normal component of the
164                 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=escript.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=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)  
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=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()  
204    
205     def __inner_PCG_v(self,v,r):       self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale
206        return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])           if self.useVPIteration:
207                  if  self.solveForFlux:
208     def __Aprod_p(self,dp):             self.__pde_k.setValue(D=self.__permeability_inv)
209        if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"          else:
210        Gdp=util.grad(dp)             self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
211        self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())          self.__pde_p.setValue(A=self.__permeability)
212        du=self.__pde_k.getSolution()           else:
213        # self.__pde_v.getOperator().saveMM("proj.mm")              D=self.__pde_k.createCoefficient("D")
214        return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)              A=self.__pde_k.createCoefficient("A")
215                D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv
216                for i in range(self.domain.getDim()):
217                   for j in range(self.domain.getDim()):
218                     A[i,i,j,j]=self.__l2
219                A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability
220                self.__pde_k.setValue(A=A, D=D)
221          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=escript.Data(), r=escript.Data(), y=escript.Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
   
    def __L2(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))  
   
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 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.3452

  ViewVC Help
Powered by ViewVC 1.1.26