/[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 2156 by gross, Mon Dec 15 05:09:02 2008 UTC revision 2251 by gross, Fri Feb 6 06:50:39 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):      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, show_details=False):      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=fixed_flux)          g=self.__g
144            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:
149               self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
150          return self.__pde_v.getSolution(verbose=show_details)          return self.__pde_v.getSolution(verbose=show_details)
151    
152      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      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            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:
215               f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
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           u2=u0*self.__pde_v.getCoefficient("q")           ATOL=self.getTolerance()
275             if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276             #########################################################################################################################
277           #           #
278           # first the reference velocity is calculated from           #   we solve:
279           #           #  
280           #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
          #  
          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)  
          u_ref=self.__pde_v.getSolution(verbose=show_details)  
          if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)  
          self.__pde_v.setValue(r=Data())  
          #  
          #   and then we calculate a reference pressure  
281           #           #
282           #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)           #   residual is
283           #           #
284           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)           #    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)
          p_ref=self.__pde_p.getSolution(verbose=self.show_details)  
          if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)  
          self.__pde_p.setValue(r=Data())  
285           #           #
286           #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
          #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref  
287           #           #
288           #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du           #    we use (g - Qp, v) to represent the residual. not that
289           #           #
290           #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
          #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref  
291           #           #
292           #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)           #   while the initial residual is
          #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)  
293           #           #
294           #  we use the (du,Qdp) to represent the resudual           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295           #  Q^*Q is a preconditioner           #  
296           #             d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297           #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>           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           Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300           norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))           self.__pde_v.setValue(r=Data())
301           ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)           # start CG
302           if not ATOL>0:           r=ArithmeticTuple(d0, v00)
303               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, norm_rhs, atol)           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           if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)           return r[1],p
305           #  
306           #   caclulate the initial residual      def __Aprod_PCG(self,dp):
          #  
          self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())  
          du=self.__pde_v.getSolution(verbose=show_details)  
          r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)  
          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)  
          util.saveVTK("d.vtu",p=dp,p_ref=p_ref)  
          return u_ref+r[1],dp  
           
     def __Aprod_PCG(self,p):  
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)
309            self.__pde_v.setValue(Y=Qp,X=Data())            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
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):
315           a=util.tensor_mult(self.__permeability,util.grad(p))           a=util.tensor_mult(self.__permeability,util.grad(p))
316           out=-util.integrate(util.inner(a,r[0]+r[1]))           f0=util.integrate(util.inner(a,r[0]))
317           return out           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    
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 264  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 294  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 324  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):
         returns div(v)  
         @rtype: equal to the type of p  
   
         @note: boundary conditions on p should be zero!  
         """  
         if self.show_details: print "apply divergence:"  
         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):  
398           """           """
399           returns inner product of element p and Bv  (overwrite)           returns inner product of element p and div(v)
           
          @type p: equal to the type of p  
          @type Bv: equal to the type of result of operator B  
          @rtype: C{float}  
400    
401           @rtype: equal to the type of p           @param p: a pressure increment
402             @param v: a residual
403             @return: inner product of element p and div(v)
404             @rtype: C{float}
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)
446             return  out
447    
448    
449             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)           out=self.__pde_u.getSolution(verbose=self.show_details)
472           return  out           return  out
473    
474        def solve_prec(self,p):       def solve_precB(self,v):
475           if self.show_details: print "apply preconditioner:"           """
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.2156  
changed lines
  Added in v.2251

  ViewVC Help
Powered by ViewVC 1.1.26