/[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 2108 by gross, Fri Nov 28 05:09:23 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 137  class DarcyFlow(object): Line 137  class DarcyFlow(object):
137                 for the permeability M{k_{ij}}                 for the permeability M{k_{ij}}
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=boundary_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")
189           g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))           #
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             #
224             Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
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             #
233             self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
234             du=self.__pde_v.getSolution(verbose=show_details)
235             r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
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"
242            Qp=util.tensor_mult(self.__permeability,util.grad(p))            Qp=util.tensor_mult(self.__permeability,util.grad(p))
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):
248           a=util.tensor_mult(self.__permeability,util.grad(p))           a=util.tensor_mult(self.__permeability,util.grad(p))
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    
257  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
# Line 257  class StokesProblemCartesian(Homogeneous Line 291  class StokesProblemCartesian(Homogeneous
291                            
292           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
293           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
294           self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
295           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
296    
297           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)

Legend:
Removed from v.2108  
changed lines
  Added in v.2156

  ViewVC Help
Powered by ViewVC 1.1.26