/[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 2208 by gross, Mon Jan 12 06:37:07 2009 UTC revision 2251 by gross, Fri Feb 6 06:50:39 2009 UTC
# Line 324  class DarcyFlow(object): Line 324  class DarcyFlow(object):
324            return self.__pde_p.getSolution(verbose=self.show_details)            return self.__pde_p.getSolution(verbose=self.show_details)
325    
326  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
327        """       """
328        solves       solves
329    
330            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
331                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 333  class StokesProblemCartesian(Homogeneous
333            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
334            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
335    
336        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
337    
338        typical usage:       typical usage:
339    
340              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
341              sp.setTolerance()              sp.setTolerance()
342              sp.initialize(...)              sp.initialize(...)
343              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
344        """       """
345        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
346           """           """
347           initialize the Stokes Problem           initialize the Stokes Problem
348    
# Line 363  class StokesProblemCartesian(Homogeneous Line 363  class StokesProblemCartesian(Homogeneous
363           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
365    
366           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
          self.__pde_proj.setReducedOrderOn()  
          self.__pde_proj.setSymmetryOn()  
          self.__pde_proj.setValue(D=1.)  
   
       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):  
367          """          """
368          assigns values to the model parameters          assigns values to the model parameters
369    
# Line 393  class StokesProblemCartesian(Homogeneous Line 388  class StokesProblemCartesian(Homogeneous
388              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
389              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
390      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
391          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
392            self.__f=f
393            self.__surface_stress=surface_stress
394          self.__stress=stress          self.__stress=stress
395    
396        def B(self,v):  #===============================================================================================================
397          """       def inner_pBv(self,p,v):
398          returns div(v)           """
399          @rtype: equal to the type of p           returns inner product of element p and div(v)
400    
401          @note: boundary conditions on p should be zero!           @param p: a pressure increment
402          """           @param v: a residual
403          if self.show_details: print "apply divergence:"           @return: inner product of element p and div(v)
         self.__pde_proj.setValue(Y=-util.div(v))  
         self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
         return self.__pde_proj.getSolution(verbose=self.show_details)  
   
       def inner_pBv(self,p,Bv):  
          """  
          returns inner product of element p and Bv  (overwrite)  
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
404           @rtype: C{float}           @rtype: C{float}
   
          @rtype: equal to the type of p  
405           """           """
406           s0=util.interpolate(p,Function(self.domain))           return util.integrate(-p*util.div(v))
          s1=util.interpolate(Bv,Function(self.domain))  
          return util.integrate(s0*s1)  
407    
408        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
409           """           """
410           returns inner product of element p0 and p1  (overwrite)           Returns inner product of p0 and p1
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
411    
412           @rtype: equal to the type of p           @param p0: a pressure
413             @param p1: a pressure
414             @return: inner product of p0 and p1
415             @rtype: C{float}
416           """           """
417           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
418           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
419           return util.integrate(s0*s1)           return util.integrate(s0*s1)
420    
421        def inner_v(self,v0,v1):       def norm_v(self,v):
422           """           """
423           returns inner product of two element v0 and v1  (overwrite)           returns the norm of v
           
          @type v0: equal to the type of v  
          @type v1: equal to the type of v  
          @rtype: C{float}  
424    
425           @rtype: equal to the type of v           @param v: a velovity
426             @return: norm of v
427             @rtype: non-negative C{float}
428           """           """
429       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))))
      gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
430    
431        def solve_A(self,u,p):       def getV(self, p, v0):
432           """           """
433           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
434    
435             @param p: a pressure
436             @param v0: a initial guess for the value v to return.
437             @return: v given as M{v= A^{-1} (f-B^*p)}
438           """           """
          if self.show_details: print "solve for velocity:"  
439           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
440             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
441           if self.__stress.isEmpty():           if self.__stress.isEmpty():
442              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain))
443           else:           else:
444              self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
445           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
446           return  out           return  out
447    
448        def solve_prec(self,p):  
449           if self.show_details: print "apply preconditioner:"           raise NotImplementedError,"no v calculation implemented."
450    
451    
452         def norm_Bv(self,v):
453            """
454            Returns Bv (overwrite).
455    
456            @rtype: equal to the type of p
457            @note: boundary conditions on p should be zero!
458            """
459            return util.sqrt(util.integrate(util.div(v)**2))
460    
461         def solve_AinvBt(self,p):
462             """
463             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
464    
465             @param p: a pressure increment
466             @return: the solution of M{Av=B^*p}
467             @note: boundary conditions on v should be zero!
468             """
469             self.__pde_u.setTolerance(self.getSubProblemTolerance())
470             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
471             out=self.__pde_u.getSolution(verbose=self.show_details)
472             return  out
473    
474         def solve_precB(self,v):
475             """
476             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
477             with accuracy L{self.getSubProblemTolerance()} (overwrite).
478    
479             @param v: velocity increment
480             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
481             @note: boundary conditions on p are zero.
482             """
483             self.__pde_prec.setValue(Y=-util.div(v))
484           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
485           self.__pde_prec.setValue(Y=p)           return self.__pde_prec.getSolution(verbose=self.show_details)
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  

Legend:
Removed from v.2208  
changed lines
  Added in v.2251

  ViewVC Help
Powered by ViewVC 1.1.26