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

revision 2123 by gross, Wed Dec 3 03:26:02 2008 UTC revision 2156 by gross, Mon Dec 15 05:09:02 2008 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 119  class DarcyFlow(object): Line 119  class DarcyFlow(object):
119             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))
120
121
122      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
123          """          """
124          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}
125          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}).
# Line 138  class DarcyFlow(object): Line 138  class DarcyFlow(object):
138          """          """
139          self.__pde_v.setTolerance(tol)          self.__pde_v.setTolerance(tol)
140          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
141          return self.__pde_v.getSolution()          return self.__pde_v.getSolution(verbose=show_details)
142
143      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
144           """           """
# Line 185  class DarcyFlow(object): Line 185  class DarcyFlow(object):
185           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
186           self.__pde_v.setTolerance(sub_rtol)           self.__pde_v.setTolerance(sub_rtol)
187           self.__pde_p.setTolerance(sub_rtol)           self.__pde_p.setTolerance(sub_rtol)
p2=p0*self.__pde_p.getCoefficient("q")
188           u2=u0*self.__pde_v.getCoefficient("q")           u2=u0*self.__pde_v.getCoefficient("q")
190           f=self.__f-util.div(u2)           # first the reference velocity is calculated from
191           self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())           #
192           dv=self.__pde_v.getSolution(verbose=show_details)           #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
193           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))           #
194             self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
195             u_ref=self.__pde_v.getSolution(verbose=show_details)
196             if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
197             self.__pde_v.setValue(r=Data())
198             #
199             #   and then we calculate a reference pressure
200             #
201             #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
202             #
203             self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
204             p_ref=self.__pde_p.getSolution(verbose=self.show_details)
205             if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
206           self.__pde_p.setValue(r=Data())           self.__pde_p.setValue(r=Data())
207           dp=self.__pde_p.getSolution(verbose=self.show_details)           #
208           norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))           #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
209           if norm_rhs<0:           #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
210               raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."           #
211           ATOL=util.sqrt(norm_rhs)*rtol +atol           #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
212             #
213             #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
214             #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
215             #
216             #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
217             #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
218             #
219             #  we use the (du,Qdp) to represent the resudual
220             #  Q^*Q is a preconditioner
221             #
222             #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
223             #
225             norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
226             ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
227           if not ATOL>0:           if not ATOL>0:
228               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, norm_rhs, atol)
229           rhs=ArithmeticTuple(g,dv)           if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
230           dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0-p2, verbose=self.verbose, initial_guess=True)           #
231           return u2+r[1],p2+dp           #   caclulate the initial residual
232             #
234             du=self.__pde_v.getSolution(verbose=show_details)
236             dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
237             util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
238             return u_ref+r[1],dp
239
240      def __Aprod_PCG(self,p):      def __Aprod_PCG(self,p):
241            if self.show_details: print "DarcyFlux: Applying operator"            if self.show_details: print "DarcyFlux: Applying operator"
243            self.__pde_v.setValue(Y=Qp,X=Data())            self.__pde_v.setValue(Y=Qp,X=Data())
244            w=self.__pde_v.getSolution(verbose=self.show_details)            w=self.__pde_v.getSolution(verbose=self.show_details)
245            return ArithmeticTuple(Qp,w)            return ArithmeticTuple(-Qp,w)
246
247      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
249           return util.integrate(util.inner(a,r[0]-r[1]))           out=-util.integrate(util.inner(a,r[0]+r[1]))
250             return out
251
252      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
253            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
254            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]))
255            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
256