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

revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC
# Line 34  __author__="Lutz Gross, l.gross@uq.edu.a Line 34  __author__="Lutz Gross, l.gross@uq.edu.a
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
38
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
# Line 48  class DarcyFlow(object): Line 48  class DarcyFlow(object):
48      @note: The problem is solved in a least squares formulation.      @note: The problem is solved in a least squares formulation.
49      """      """
50
51      def __init__(self, domain):      def __init__(self, domain,useReduced=False):
52          """          """
53          initializes the Darcy flux problem          initializes the Darcy flux problem
54          @param domain: domain of the problem          @param domain: domain of the problem
# Line 56  class DarcyFlow(object): Line 56  class DarcyFlow(object):
56          """          """
57          self.domain=domain          self.domain=domain
58          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
59          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          if useReduced: self.__pde_v.setReducedOrderOn()
60          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
61            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
62          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
63          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
64            if useReduced: self.__pde_p.setReducedOrderOn()
65          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67            self.__ATOL= None
68
69      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):
70          """          """
# Line 119  class DarcyFlow(object): Line 122  class DarcyFlow(object):
122             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
123
124
125      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):      def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
126          """          """
127          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
128          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
129          Note that C{g} and C{f} are used, L{setValue}.          Note that C{g} and C{f} are used, see L{setValue}.
130
131          @param p: pressure.          @param p: pressure.
132          @type p: scalar value on the domain (e.g. L{Data}).          @type p: scalar value on the domain (e.g. L{Data}).
133          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134          @type fixed_flux: vector values on the domain (e.g. L{Data}).          @type fixed_flux: vector values on the domain (e.g. L{Data}).
135          @param tol: relative tolerance to be used.          @param tol: relative tolerance to be used.
136          @type tol: positive float.          @type tol: positive C{float}.
137          @return: flux          @return: flux
138          @rtype: L{Data}          @rtype: L{Data}
139          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
140                 for the permeability M{k_{ij}}                 for the permeability M{k_{ij}}
141          """          """
142          self.__pde_v.setTolerance(tol)          self.__pde_v.setTolerance(tol)
143          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)          g=self.__g
144          return self.__pde_v.getSolution()          f=self.__f
145            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
146            if p == None:
147               self.__pde_v.setValue(Y=g)
148            else:
150            return self.__pde_v.getSolution(verbose=show_details)
151
152        def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
153            """
154            returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
155            on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).
156            Note that C{g} is used, see L{setValue}.
157
158            @param v: flux.
159            @type v: vector-valued on the domain (e.g. L{Data}).
160            @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
161            @type fixed_pressure: vector values on the domain (e.g. L{Data}).
162            @param tol: relative tolerance to be used.
163            @type tol: positive C{float}.
164            @return: pressure
165            @rtype: L{Data}
166            @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}
167                   for the permeability M{k_{ij}}
168            """
169            self.__pde_v.setTolerance(tol)
170            g=self.__g
171            self.__pde_p.setValue(r=fixed_pressure)
172            if v == None:
173               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
174            else:
175               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
176            return self.__pde_p.getSolution(verbose=show_details)
177
178        def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
179            """
180            set the tolerance C{ATOL} used to terminate the solution process. It is used
181
182            M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}
183
184      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):          where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.
185
186            The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual
187
188            M{r=Q^*(g-Qp-v)}
189
190            the condition
191
192            M{<(Q^*Q)^{-1} r,r> <= ATOL}
193
194            holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
195
196            @param atol: absolute tolerance for the pressure
197            @type atol: non-negative C{float}
198            @param rtol: relative tolerance for the pressure
199            @type rtol: non-negative C{float}
200            @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the
201            L{getPressure} method or use  the value from a previous time step.
202            @type p_ref: scalar value on the domain (e.g. L{Data}).
203            @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the
204            L{getFlux} method or use  the value from a previous time step.
205            @type v_ref: vector-valued on the domain (e.g. L{Data}).
206            @return: used absolute tolerance.
207            @rtype: positive C{float}
208            """
209            g=self.__g
210            if not v_ref == None:
211               f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212            else:
213               f1=util.integrate(util.length(util.interpolate(g))**2)
214            if not p_ref == None:
216            else:
217               f2=0
218            self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219            if self.__ATOL<=0:
220               raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221            return self.__ATOL
222
223        def getTolerance(self):
224            """
225            returns the current tolerance.
226
227            @return: used absolute tolerance.
228            @rtype: positive C{float}
229            """
230            if self.__ATOL==None:
231               raise ValueError,"no tolerance is defined."
232            return self.__ATOL
233
234        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235           """           """
236           solves the problem.           solves the problem.
237
238           The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where           The iteration is terminated if the residual norm is less then self.getTolerance().
C{|q|} denotes the norm of the right hand side (see escript user's guide for details).
239
240           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.           @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
241           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
242           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.           @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
243           @type p0: scalar value on the domain (e.g. L{Data}).           @type p0: scalar value on the domain (e.g. L{Data}).
@param atol: absolute tolerance for the pressure
@type atol: non-negative C{float}
@param rtol: relative tolerance for the pressure
@type rtol: non-negative C{float}
244           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245           @type sub_rtol: positive-negative C{float}           @type sub_rtol: positive-negative C{float}
246           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
# Line 185  class DarcyFlow(object): Line 271  class DarcyFlow(object):
271           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
272           self.__pde_v.setTolerance(sub_rtol)           self.__pde_v.setTolerance(sub_rtol)
273           self.__pde_p.setTolerance(sub_rtol)           self.__pde_p.setTolerance(sub_rtol)
274           p2=p0*self.__pde_p.getCoefficient("q")           ATOL=self.getTolerance()
275           u2=u0*self.__pde_v.getCoefficient("q")           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
277           f=self.__f-util.div(u2)           #
278           self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())           #   we solve:
279           dv=self.__pde_v.getSolution(verbose=show_details)           #
280           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
281           self.__pde_p.setValue(r=Data())           #
282           dp=self.__pde_p.getSolution(verbose=self.show_details)           #   residual is
283           norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))           #
284           if norm_rhs<0:           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)
285               raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."           #
286           ATOL=util.sqrt(norm_rhs)*rtol +atol           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
287           if not ATOL>0:           #
288               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)           #    we use (g - Qp, v) to represent the residual. not that
289           rhs=ArithmeticTuple(g,dv)           #
290           dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0, verbose=self.verbose, initial_guess=False)           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291           return u2+r[1],p2+dp           #
292                     #   while the initial residual is
293      def __Aprod_PCG(self,p):           #
294             #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295             #
297             self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
298             v00=self.__pde_v.getSolution(verbose=show_details)
299             if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300             self.__pde_v.setValue(r=Data())
301             # start CG
302             r=ArithmeticTuple(d0, v00)
303             p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
304             return r[1],p
305
306        def __Aprod_PCG(self,dp):
307            if self.show_details: print "DarcyFlux: Applying operator"            if self.show_details: print "DarcyFlux: Applying operator"
308            Qp=util.tensor_mult(self.__permeability,util.grad(p))            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
310            w=self.__pde_v.getSolution(verbose=self.show_details)            self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311            return ArithmeticTuple(Qp,w)            du=self.__pde_v.getSolution(verbose=self.show_details)
312              return ArithmeticTuple(mQdp,du)
313
314      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
316           return util.integrate(util.inner(a,r[0]-r[1]))           f0=util.integrate(util.inner(a,r[0]))
317             f1=util.integrate(util.inner(a,r[1]))
318             # print "__inner_PCG:",f0,f1,"->",f0-f1
319             return f0-f1
320
321      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
322            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
323            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]))            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
324            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
325