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

revision 2344 by jfenwick, Mon Mar 30 02:13:58 2009 UTC revision 2445 by gross, Fri May 29 03:23:25 2009 UTC
# 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,useReduced=False):      def __init__(self, domain, weight=None, 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
55          @type domain: L{Domain}          @type domain: L{Domain}
56          """          """
57          self.domain=domain          self.domain=domain
58          self.__l=util.longestEdge(self.domain)**2          if weight == None:
59               s=self.domain.getSize()
60               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61            else:
62               self.__l=weight
63          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
64          if useReduced: self.__pde_v.setReducedOrderOn()          if useReduced: self.__pde_v.setReducedOrderOn()
65          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
66          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.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
69          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
70          if useReduced: self.__pde_p.setReducedOrderOn()          if useReduced: self.__pde_p.setReducedOrderOn()
# Line 230  class DarcyFlow(object): Line 235  class DarcyFlow(object):
235           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
236           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237           """           """
238           self.verbose=verbose or True           self.verbose=verbose
239           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
240           rtol=self.getTolerance()           rtol=self.getTolerance()
241           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
# Line 263  class DarcyFlow(object): Line 268  class DarcyFlow(object):
268                 if norm_r == None or norm_r>ATOL:                 if norm_r == None or norm_r>ATOL:
269                     if num_corrections>max_num_corrections:                     if num_corrections>max_num_corrections:
270                           raise ValueError,"maximum number of correction steps reached."                           raise ValueError,"maximum number of correction steps reached."
271                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272                     num_corrections+=1                     num_corrections+=1
273                 else:                 else:
274                     converged=True                     converged=True
275           return v,p           return v,p
#
#
#               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
#               #===========================================================================
#               norm_r_hat=self.__L2(r_hat)
#               norm_v=self.__L2(v)
#               norm_g=self.__L2(g)
#               norm_gv=self.__L2(g-v)
#               norm_Qp=self.__L2(Qp)
#               norm_gQp=self.__L2(g-Qp)
#               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
#               fac=min(norm_v,norm_Qp,norm_gv)
#               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
#               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
#               if r!=None:
#                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
#                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
#                   self.setSubProblemTolerance(sub_tol)
#                   print "subtol_new=",self.getSubProblemTolerance()
#               print "norm_v = ",norm_v
#               print "norm_gv = ",norm_gv
#               print "norm_Qp = ",norm_Qp
#               print "norm_gQp = ",norm_gQp
#               print "norm_g = ",norm_g
#               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
#               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
#               if fac == 0:
#                   if self.verbose: print "DarcyFlux: trivial case!"
#                   return v,p
#               #===============================================================================
#               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
#               # norm_Qp=self.__L2(Qp)
#               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
#               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
#
#               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
#               if self.verbose:
#                   print "DarcyFlux: residual = %e"%norm_r_hat
#                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
#               if norm_r_hat <= ATOL:
#                   print "DarcyFlux: iteration finalized."
#                   converged=True
#               else:
#                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)
#                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)
#                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
#               print "norm_r =",norm_r
#         return v,p
276      def __L2(self,v):      def __L2(self,v):
277           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278
# Line 327  class DarcyFlow(object): Line 284  class DarcyFlow(object):
284            if self.show_details: print "DarcyFlux: Applying operator"            if self.show_details: print "DarcyFlux: Applying operator"
285            Qdp=self.__Q(dp)            Qdp=self.__Q(dp)
286            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287            du=self.__pde_v.getSolution(verbose=self.show_details)            du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288              # self.__pde_v.getOperator().saveMM("proj.mm")
289            return Qdp+du            return Qdp+du
290      def __inner_GMRES(self,r,s):      def __inner_GMRES(self,r,s):
291           return util.integrate(util.inner(r,s))           return util.integrate(util.inner(r,s))
# Line 339  class DarcyFlow(object): Line 297  class DarcyFlow(object):
297            self.__pde_p.setTolerance(self.getSubProblemTolerance())            self.__pde_p.setTolerance(self.getSubProblemTolerance())
298            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
299            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
301              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302
303      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304          """          """
# Line 366  class DarcyFlow(object): Line 325  class DarcyFlow(object):
325             self.__pde_v.setValue(Y=g)             self.__pde_v.setValue(Y=g)
326          else:          else:
327             self.__pde_v.setValue(Y=g-self.__Q(p))             self.__pde_v.setValue(Y=g-self.__Q(p))
328          return self.__pde_v.getSolution(verbose=show_details)          return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329
331       """       """
# Line 401  class StokesProblemCartesian(Homogeneous Line 360  class StokesProblemCartesian(Homogeneous
360           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
361           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
362           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364
365           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
366           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
367           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
369
370             self.__pde_proj=LinearPDE(domain)
371             self.__pde_proj.setReducedOrderOn()
372         self.__pde_proj.setValue(D=1)
373             self.__pde_proj.setSymmetryOn()
374
375
377          """          """
378          assigns values to the model parameters          assigns values to the model parameters
# Line 438  class StokesProblemCartesian(Homogeneous Line 403  class StokesProblemCartesian(Homogeneous
403          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
404          self.__stress=stress          self.__stress=stress
405
406       def inner_pBv(self,p,v):       def Bv(self,v):
407           """           """
408           returns inner product of element p and div(v)           returns inner product of element p and div(v)
409
# Line 447  class StokesProblemCartesian(Homogeneous Line 412  class StokesProblemCartesian(Homogeneous
412           @return: inner product of element p and div(v)           @return: inner product of element p and div(v)
413           @rtype: C{float}           @rtype: C{float}
414           """           """
415           return util.integrate(-p*util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
416             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
417             return self.__pde_proj.getSolution()
418
419         def inner_pBv(self,p,Bv):
420             """
421             returns inner product of element p and Bv=-div(v)
422
423             @param p: a pressure increment
424             @param v: a residual
425             @return: inner product of element p and Bv=-div(v)
426             @rtype: C{float}
427             """
428             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
429
430       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
431           """           """
# Line 489  class StokesProblemCartesian(Homogeneous Line 467  class StokesProblemCartesian(Homogeneous
467           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
468           return  out           return  out
469
470         def norm_Bv(self,Bv):
raise NotImplementedError,"no v calculation implemented."

def norm_Bv(self,v):
471          """          """
472          Returns Bv (overwrite).          Returns Bv (overwrite).
473
474          @rtype: equal to the type of p          @rtype: equal to the type of p
475          @note: boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
476          """          """
477          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
478
479       def solve_AinvBt(self,p):       def solve_AinvBt(self,p):
480           """           """
# Line 515  class StokesProblemCartesian(Homogeneous Line 489  class StokesProblemCartesian(Homogeneous
489           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
490           return  out           return  out
491
492       def solve_precB(self,v):       def solve_prec(self,Bv):
493           """           """
494           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
495           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy L{self.getSubProblemTolerance()}
496
497           @param v: velocity increment           @param v: velocity increment
498           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
499           @note: boundary conditions on p are zero.           @note: boundary conditions on p are zero.
500           """           """
501           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
502           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
503           return self.__pde_prec.getSolution(verbose=self.show_details)           return self.__pde_prec.getSolution(verbose=self.show_details)
504
505    # rename solve_prec and change argument v to Bv
506    # chnage the argument of inner_pBv to v->Bv