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

revision 3070 by gross, Tue Mar 2 07:54:11 2010 UTC revision 3071 by gross, Wed Jul 21 05:37:30 2010 UTC
# Line 1  Line 1
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
# Line 37  from linearPDEs import LinearPDE, Linear Line 38  from linearPDEs import LinearPDE, Linear
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39
40  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      solves the problem     solves the problem
43
44      *u_i+k_{ij}*p_{,j} = g_i*     *u_i+k_{ij}*p_{,j} = g_i*
45      *u_{i,i} = f*     *u_{i,i} = f*
46
47      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48
49      :note: The problem is solved in a least squares formulation.     :note: The problem is solved in a least squares formulation.
50      """     """
51
52      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):
53          """        """
54          initializes the Darcy flux problem        initializes the Darcy flux problem
55          :param domain: domain of the problem        :param domain: domain of the problem
56          :type domain: `Domain`        :type domain: `Domain`
57          :param weight: the weighting factor for the incempressiblity equation *u_{i,i} = f* in the variational principle.        :param useReduced: uses reduced oreder on flux and pressure
58          If not present an apprppriate weight is chosen.        :type useReduced: ``bool``
59          :type  weight: positive sacalar        :param adaptSubTolerance: switches on automatic subtolerance selection
60      :param useReduced: uses reduced oreder on flux and pressure        :type adaptSubTolerance: ``bool``
61      :type useReduced: ``bool``        :param solveForFlux: if True the solver solves for the flux (do not use!)
62      :param adaptSubTolerance: switches on automatic subtolerance selection        :type solveForFlux: ``bool``
63      :type adaptSubTolerance: ``bool``          """
64          """        self.domain=domain
65          self.domain=domain        self.solveForFlux=solveForFlux
66          if weight == None:        self.useReduced=useReduced
68             self.__update_weight=True        self.verbose=False
69          else:
70             self.__update_weight=False        self.__pde_k=LinearPDESystem(domain)
71             self.__l=weight        self.__pde_k.setSymmetryOn()
72          self.__pde_v=LinearPDESystem(domain)        if self.useReduced: self.__pde_k.setReducedOrderOn()
73          if useReduced: self.__pde_v.setReducedOrderOn()
74          self.__pde_v.setSymmetryOn()        self.__pde_p=LinearSinglePDE(domain)
75          self.__pde_p=LinearSinglePDE(domain)        self.__pde_p.setSymmetryOn()
76          self.__pde_p.setSymmetryOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
77          if useReduced: self.__pde_p.setReducedOrderOn()
78          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting
79          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        # self.__pde_l.setSymmetryOn()
80          self.setTolerance()        # if self.useReduced: self.__pde_l.setReducedOrderOn()
81          self.setAbsoluteTolerance()        self.setTolerance()
83      self.verbose=False        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84      def getSolverOptionsFlux(self):        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85      """
86      Returns the solver options used to solve the flux problems     def getSolverOptionsFlux(self):
87              """
88      *(I+D^* weight D)u=F*        Returns the solver options used to solve the flux problems
89
90      :return: `SolverOptions`        *K^{-1} u=F*
91      """
92      return self.__pde_v.getSolverOptions()        :return: `SolverOptions`
93          """
94      def setSolverOptionsFlux(self, options=None):        return self.__pde_k.getSolverOptions()
95      """
96      Sets the solver options used to solve the flux problems     def setSolverOptionsFlux(self, options=None):
97              """
98      *(I+D^*  weight  D)u=F*        Sets the solver options used to solve the flux problems
99
100      If ``options`` is not present, the options are reset to default        *K^{-1}u=F*
101      :param options: `SolverOptions`
102      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.        If ``options`` is not present, the options are reset to default
103      """
104      return self.__pde_v.setSolverOptions(options)        :param options: `SolverOptions`
105      def getSolverOptionsPressure(self):        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
106      """        """
107      Returns the solver options used to solve the pressure problems        return self.__pde_v.setSolverOptions(options)
108
109      *(Q^* S Q)p=Q^*G*     def getSolverOptionsPressure(self):
110              """
111           as a preconditioner where S=weight*k^2/dx^2        Returns the solver options used to solve the pressure problems
112
113      :return: `SolverOptions`        *(Q^* K Q)p=-Q^*G*
114      """
115      return self.__pde_p.getSolverOptions()        :return: `SolverOptions`
116      def setSolverOptionsPressure(self, options=None):        """
117      """        return self.__pde_p.getSolverOptions()
118      Sets the solver options used to solve the pressure problems
119           def setSolverOptionsPressure(self, options=None):
120      *(Q^* S Q)p=Q^*G*        """
121              Sets the solver options used to solve the pressure problems
122          was a preconditioner here S=weight*k^2/dx^2
123              *(Q^* K Q)p=-Q^*G*
124      If ``options`` is not present, the options are reset to default
125          If ``options`` is not present, the options are reset to default
126      :param options: `SolverOptions`
127      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.        :param options: `SolverOptions`
128      """        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
129      return self.__pde_p.setSolverOptions(options)        """
130          return self.__pde_p.setSolverOptions(options)
131      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
132          """     def getSolverOptionsWeighting(self):
133          assigns values to model parameters        """
134          Returns the solver options used to solve the pressure problems
135          :param f: volumetic sources/sinks
136          :type f: scalar value on the domain (e.g. `Data`)        *(D K D^*)p=-D F*
137          :param g: flux sources/sinks
138          :type g: vector values on the domain (e.g. `Data`)        :return: `SolverOptions`
139          :param location_of_fixed_pressure: mask for locations where pressure is fixed        """
140          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)        return self.__pde_l.getSolverOptions()
141          :param location_of_fixed_flux:  mask for locations where flux is fixed.
142          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)     def setSolverOptionsWeighting(self, options=None):
143          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with        """
144                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with        Sets the solver options used to solve the pressure problems
145                               ``v`` on the main diagonal is used.
146          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)        *(D K D^*)p=-D F*
147
148          :note: the values of parameters which are not set by calling ``setValue`` are not altered.        If ``options`` is not present, the options are reset to default
149          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
150                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal        :param options: `SolverOptions`
151                 is along the *x_i* axis.        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
152          """        """
153          if f !=None:        return self.__pde_l.setSolverOptions(options)
154             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
155             if f.isEmpty():
156                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
157             else:        """
158                 if f.getRank()>0: raise ValueError,"illegal rank of f."        assigns values to model parameters
159             self.__f=f
160          if g !=None:        :param f: volumetic sources/sinks
161             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))        :type f: scalar value on the domain (e.g. `Data`)
162             if g.isEmpty():        :param g: flux sources/sinks
163               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        :type g: vector values on the domain (e.g. `Data`)
164             else:        :param location_of_fixed_pressure: mask for locations where pressure is fixed
165               if not g.getShape()==(self.domain.getDim(),):        :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
166                 raise ValueError,"illegal shape of g"        :param location_of_fixed_flux:  mask for locations where flux is fixed.
167             self.__g=g        :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
168          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
169          if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)        ``s`` on the main diagonal is used.
170          if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)        :type permeability: scalar or tensor values on the domain (e.g. `Data`)
171          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
172          if permeability!=None:        :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
173             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))        or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
174             if perm.getRank()==0:        is along the *x_i* axis.
175                 perm_inv=1./perm*util.kronecker(self.domain.getDim())        """
176                 perm=perm*util.kronecker(self.domain.getDim())        if f !=None:
177             elif perm.getRank()==1:       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178                 perm2=perm       if f.isEmpty():
179                 perm=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A"))          f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180                 perm_inv=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A"))       else:
181                 for i in range(self.domain.getDim()):          if f.getRank()>0: raise ValueError,"illegal rank of f."
182                        p=perm2[i]          self.__f=f
183                        perm[i,i]=p        if g !=None:
184                        perm_inv[i,i]=1/p       g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185             elif perm.getRank()==2:       if g.isEmpty():
186                 perm_inv=util.inverse(perm)          g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187             else:       else:
188                raise ValueError,"illegal rank of permeability."          if not g.getShape()==(self.domain.getDim(),):
189             self.__permeability=perm             raise ValueError,"illegal shape of g"
190             self.__permeability_inv=perm_inv          self.__g=g
191             if self.__update_weight: self.__l = util.longestEdge(self.domain)**2/3.14**2*util.length(self.__permeability_inv)        if location_of_fixed_pressure!=None:
192             self.__pde_v.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))             self.__pde_p.setValue(q=location_of_fixed_pressure)
193             s=self.__pde_p.getFunctionSpaceForCoefficient("A").getSize()             #self.__pde_l.setValue(q=location_of_fixed_pressure)
194             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)*self.__l/s**2)        if location_of_fixed_flux!=None:
195               self.__pde_k.setValue(q=location_of_fixed_flux)
196      def getFlux(self,p=None, fixed_flux=Data()):
197          """        if permeability!=None:
198          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
199          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).       if perm.getRank()==0:
200          Note that ``g`` and ``f`` are used, see `setValue`.          perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
201            perm=perm*util.kronecker(self.domain.getDim())
202          :param p: pressure.       elif perm.getRank()==2:
203          :type p: scalar value on the domain (e.g. `Data`).          perm_inv=util.inverse(perm)
204          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.       else:
205          :type fixed_flux: vector values on the domain (e.g. `Data`).          raise ValueError,"illegal rank of permeability."
206          :return: flux
207          :rtype: `Data`       self.__permeability=perm
208          :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}*       self.__permeability_inv=perm_inv
209                 for the permeability *k_{ij}*       self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
210          """       #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
211      self.setSubProblemTolerance()       if  self.solveForFlux:
212          self.__pde_v.setValue(X=util.sqrt(self.__l)*self.__f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_k.setValue(D=self.__permeability_inv)
213          g2=util.tensor_mult(self.__permeability_inv,self.__g)       else:
214          if p == None:          self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
215             self.__pde_v.setValue(Y=g2)       self.__pde_p.setValue(A=self.__permeability)
216          else:       #self.__pde_l.setValue(D=1/self.__l)
218          return self.__pde_v.getSolution()
219       def __applWeight(self, v, f=None):
220      def __Aprod(self,dp):        # solves L p = f-Dv with p = 0
221            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"        if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
222            Gdp=util.grad(dp)        if f == None:
223            self.__pde_v.setValue(Y=-Gdp,X=Data(), r=Data())       return -util.div(v)*self.__l
224            du=self.__pde_v.getSolution()        else:
225            # self.__pde_v.getOperator().saveMM("proj.mm")       return (f-util.div(v))*self.__l
226            return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)        # if f == None:
227          #      self.__pde_l.setValue(Y=-util.div(v))
228      def __Msolve_PCG(self,r):        # else:
229        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        #      return (f-util.div(v))/self.__l
230            self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data())        # return self.__pde_l.getSolution()
231            # self.__pde_p.getOperator().saveMM("prec.mm")
232            return self.__pde_p.getSolution()     def __getPressure(self, v, p0, g=None):
233          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
234      def __inner_PCG(self,p,r):        if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
235           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))        if g == None:
236         self.__pde_p.setValue(X=-v, r=p0)
237      def __L2(self,v):        else:
238           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))       self.__pde_p.setValue(X=g-v, r=p0)
239          p=self.__pde_p.getSolution()
240      def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):        return p
241           """
242           solves the problem.     def __Aprod_v(self,dv):
243          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv
244           The iteration is terminated if the residual norm is less then self.getTolerance().        dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
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.        b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247           :type u0: vector value on the domain (e.g. `Data`).        return ArithmeticTuple(a,b,-dp)
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.
249           :type p0: scalar value on the domain (e.g. `Data`).     def __Msolve_PCG_v(self,r):
250           :param verbose: if set some information on iteration progress are printed        # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251           :type verbose: ``bool``        if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252           :return: flux and pressure        self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())
253           :rtype: ``tuple`` of `Data`.        # self.__pde_p.getOperator().saveMM("prec.mm")
254          return self.__pde_k.getSolution()
255           :note: The problem is solved as a least squares form
256       def __inner_PCG_v(self,v,r):
257           *(K^[-1]+D^* weight D)u+G p=D^* sqrt(weight) f + K^[-1]g*        return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
258           *G^*u+G^* K Gp=G^*g*
259       def __Aprod_p(self,dp):
260           where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.        if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261           We eliminate the flux form the problem by setting        Gdp=util.grad(dp)
262          self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())
263           *u=(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g* -G p ) with u=u0 on location_of_fixed_flux        du=self.__pde_k.getSolution()
264          # self.__pde_v.getOperator().saveMM("proj.mm")
265           form the first equation. Inserted into the second equation we get        return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
266
267           *G^*(K-(K^[-1]+D^* weight D)^{-1})Gp= G^*(g-(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g))* with p=p0  on location_of_fixed_pressure     def __getFlux(self,p, v0, f=None, g=None):
268          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
269           which is solved using the PCG method (precondition is *G^*K^2*weight*l^2/dx^2 G*. In each iteration step        if f!=None:
270           PDEs with operator *K^[-1]+D^* weight D* and with *G^*K^2*weight*l^2/dx^2 G* needs to be solved using a sub iteration scheme.       self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
271           """        self.__pde_k.setValue(r=v0)
272           self.verbose=verbose        g2=util.tensor_mult(self.__permeability_inv,g)
273           rtol=self.getTolerance()        if p == None:
274           atol=self.getAbsoluteTolerance()       self.__pde_k.setValue(Y=g2)
275       self.setSubProblemTolerance()        else:
277           converged=False        return self.__pde_k.getSolution()
278           p=p0
279           norm_r=None        #v=self.__getFlux(p, u0, f=self.__f, g=g2)
280           v=self.getFlux(p, fixed_flux=u0)     def __Msolve_PCG_p(self,r):
281           while not converged:        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282                 Gp=util.grad(p)        self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())
283                 KGp=util.tensor_mult(self.__permeability,Gp)        # self.__pde_p.getOperator().saveMM("prec.mm")
284                 norm_v=util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5        return self.__pde_p.getSolution()
285                 norm_Gp=util.integrate(util.inner(Gp,KGp))**0.5
286                 if self.verbose:     def __inner_PCG_p(self,p,r):
287                     print "DarcyFlux: L2: g-v-K*grad(p) = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-KGp)**2)**(0.5),)         return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288                     print "DarcyFlux: L2: f-div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
289                     print "DarcyFlux: K^{-1}-norm of v = %e."%norm_v     def __L2(self,v):
290                     print "DarcyFlux: K-norm of grad(p) = %e."%norm_Gp        return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291                 if norm_v == 0.:
292                    if norm_Gp == 0.:     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293                       return v,p        """
294                    else:        solves the problem.
295                      fac=norm_Gp
296                 else:        The iteration is terminated if the residual norm is less then self.getTolerance().
297                    if norm_Gp == 0.:
298                      fac=norm_v        :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.
299                    else:        :type u0: vector value on the domain (e.g. `Data`).
300                      fac=2./(1./norm_v+1./norm_Gp)        :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.
301                 #fac=max(norm_Gp, norm_v)        :type p0: scalar value on the domain (e.g. `Data`).
302                 ATOL=atol+rtol*fac        :param verbose: if set some information on iteration progress are printed
303                 if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e (fac= %e)."%(ATOL,fac)        :type verbose: ``bool``
304                 if norm_r == None or norm_r>ATOL:        :return: flux and pressure
305                     if num_corrections>max_num_corrections:        :rtype: ``tuple`` of `Data`.
306                           raise ValueError,"maximum number of correction steps reached."
307                     # initial residual is r=G^*(self.__g-KGp - v)        :note: The problem is solved as a least squares form
308                     p,r, norm_r=PCG(ArithmeticTuple(self.__g-KGp,v),        *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
309                                     self.__Aprod,        *G^*u+G^* K Gp=G^*g*
310                                     p,
311                                     self.__Msolve_PCG,        where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
312                                     self.__inner_PCG,        """
313                                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)        self.verbose=verbose
314                     if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r        rtol=self.getTolerance()
315                     v=r[1]        atol=self.getAbsoluteTolerance()
316                     num_corrections+=1        self.setSubProblemTolerance()
317                 else:        num_corrections=0
318                     if self.verbose: print "DarcyFlux: stopping criterium reached."        converged=False
319                     converged=True        norm_r=None
320           return v,p
321          # Eliminate the hydrostatic pressure:
322                    if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
323          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))
324      def setTolerance(self,rtol=1e-4):        p0=self.__pde_p.getSolution()
325          """        g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
326          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5
327
328          *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}*        p=p0*0
329          if self.solveForFlux:
330          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).       v=u0.copy()
331          else:
332          :param rtol: relative tolerance for the pressure       v=self.__getFlux(p, u0, f=self.__f, g=g2)
333          :type rtol: non-negative ``float``
334          """        while not converged and norm_g2 > 0:
335          if rtol<0:       Gp=util.grad(p)
336              raise ValueError,"Relative tolerance needs to be non-negative."       KGp=util.tensor_mult(self.__permeability,Gp)
337          self.__rtol=rtol       if self.verbose:
338      def getTolerance(self):          def_p=g2-(v+KGp)
339          """          def_v=self.__f-util.div(v)
340          returns the relative tolerance          print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
341            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
342          :return: current relative tolerance          print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
343          :rtype: ``float``          print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
344          """          print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
345          return self.__rtol       ATOL=atol+rtol*norm_g2
346         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
347      def setAbsoluteTolerance(self,atol=0.):       if norm_r == None or norm_r>ATOL:
348          """          if num_corrections>max_num_corrections:
349          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if             raise ValueError,"maximum number of correction steps reached."
350
351           *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}*          if self.solveForFlux:
352               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
353          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.             v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
354                   self.__Aprod_v,
355          :param atol: absolute tolerance for the pressure                 v,
356          :type atol: non-negative ``float``                 self.__Msolve_PCG_v,
357          """                 self.__inner_PCG_v,
358          if atol<0:                 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
359              raise ValueError,"Absolute tolerance needs to be non-negative."             p=r[2]
360          self.__atol=atol          else:
361      def getAbsoluteTolerance(self):             # initial residual is r=G^*(g2-KGp - v)
362         """             p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
363         returns the absolute tolerance                   self.__Aprod_p,
364                           p,
365         :return: current absolute tolerance                   self.__Msolve_PCG_p,
366         :rtype: ``float``                   self.__inner_PCG_p,
367         """                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
368         return self.__atol             v=r[1]
369      def getSubProblemTolerance(self):          if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
370      """          num_corrections+=1
371      Returns a suitable subtolerance       else:
372      @type: ``float``          if self.verbose: print "DarcyFlux: stopping criterium reached."
373      """          converged=True
374      return max(util.EPSILON**(0.75),self.getTolerance()**2)        return v,p+p0
375       def setTolerance(self,rtol=1e-4):
376      def setSubProblemTolerance(self):        """
377           """        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
378           Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
379           """        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
381           sub_tol=self.getSubProblemTolerance()        where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
382               self.getSolverOptionsFlux().setTolerance(sub_tol)
383           self.getSolverOptionsFlux().setAbsoluteTolerance(0.)        :param rtol: relative tolerance for the pressure
384           self.getSolverOptionsPressure().setTolerance(sub_tol)        :type rtol: non-negative ``float``
385           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)        """
386           if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol        if rtol<0:
387         raise ValueError,"Relative tolerance needs to be non-negative."
388          self.__rtol=rtol
389       def getTolerance(self):
390          """
391          returns the relative tolerance
392          :return: current relative tolerance
393          :rtype: ``float``
394          """
395          return self.__rtol
396
397       def setAbsoluteTolerance(self,atol=0.):
398          """
399          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
400
401          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
402
403
404          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
405
406          :param atol: absolute tolerance for the pressure
407          :type atol: non-negative ``float``
408          """
409          if atol<0:
410         raise ValueError,"Absolute tolerance needs to be non-negative."
411          self.__atol=atol
412       def getAbsoluteTolerance(self):
413          """
414          returns the absolute tolerance
415          :return: current absolute tolerance
416          :rtype: ``float``
417          """
418          return self.__atol
419       def getSubProblemTolerance(self):
420          """
421          Returns a suitable subtolerance
422          :type: ``float``
423          """
424          return max(util.EPSILON**(0.75),self.getTolerance()**2)
425
426       def setSubProblemTolerance(self):
427          """
428          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
429          """
431         sub_tol=self.getSubProblemTolerance()
432         self.getSolverOptionsFlux().setTolerance(sub_tol)
433         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
434         self.getSolverOptionsPressure().setTolerance(sub_tol)
435         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
436         self.getSolverOptionsWeighting().setTolerance(sub_tol)
437         self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
438         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
439
440
441  class DarcyFlowOld(object):  class DarcyFlowOld(object):

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