/[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 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 21  __url__="https://launchpad.net/escript-f Line 21  __url__="https://launchpad.net/escript-f
21  """  """
22  Some models for flow  Some models for flow
23    
24  @var __author__: name of author  :var __author__: name of author
25  @var __copyright__: copyrights  :var __copyright__: copyrights
26  @var __license__: licence agreement  :var __license__: licence agreement
27  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
28  @var __version__: version  :var __version__: version
29  @var __date__: date of the version  :var __date__: date of the version
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 40  class DarcyFlow(object): Line 40  class DarcyFlow(object):
40      """      """
41      solves the problem      solves the problem
42    
43      M{u_i+k_{ij}*p_{,j} = g_i}      *u_i+k_{ij}*p_{,j} = g_i*
44      M{u_{i,i} = f}      *u_{i,i} = f*
45    
46      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
47    
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, weight=None, useReduced=False, adaptSubTolerance=True):      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
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: `Domain`
56      @param useReduced: uses reduced oreder on flux and pressure      :param useReduced: uses reduced oreder on flux and pressure
57      @type useReduced: C{bool}      :type useReduced: ``bool``
58      @param adaptSubTolerance: switches on automatic subtolerance selection      :param adaptSubTolerance: switches on automatic subtolerance selection
59      @type adaptSubTolerance: C{bool}          :type adaptSubTolerance: ``bool``  
60          """          """
61          self.domain=domain          self.domain=domain
62          if weight == None:          if weight == None:
63             s=self.domain.getSize()             s=self.domain.getSize()
64             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l=(3.*util.longestEdge(self.domain))**2
66             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
67          else:          else:
68             self.__l=weight             self.__l=weight
69          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
# Line 83  class DarcyFlow(object): Line 83  class DarcyFlow(object):
83      """      """
84      Returns the solver options used to solve the flux problems      Returns the solver options used to solve the flux problems
85            
86      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
87            
88      @return: L{SolverOptions}      :return: `SolverOptions`
89      """      """
90      return self.__pde_v.getSolverOptions()      return self.__pde_v.getSolverOptions()
91      def setSolverOptionsFlux(self, options=None):      def setSolverOptionsFlux(self, options=None):
92      """      """
93      Sets the solver options used to solve the flux problems      Sets the solver options used to solve the flux problems
94            
95      M{(I+D^*D)u=F}      *(I+D^*D)u=F*
96            
97      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
98      @param options: L{SolverOptions}      :param options: `SolverOptions`
99      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
100      """      """
101      return self.__pde_v.setSolverOptions(options)      return self.__pde_v.setSolverOptions(options)
102      def getSolverOptionsPressure(self):      def getSolverOptionsPressure(self):
103      """      """
104      Returns the solver options used to solve the pressure problems      Returns the solver options used to solve the pressure problems
105            
106      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
107            
108      @return: L{SolverOptions}      :return: `SolverOptions`
109      """      """
110      return self.__pde_p.getSolverOptions()      return self.__pde_p.getSolverOptions()
111      def setSolverOptionsPressure(self, options=None):      def setSolverOptionsPressure(self, options=None):
112      """      """
113      Sets the solver options used to solve the pressure problems      Sets the solver options used to solve the pressure problems
114            
115      M{(Q^*Q)p=Q^*G}      *(Q^*Q)p=Q^*G*
116            
117      If C{options} is not present, the options are reset to default      If ``options`` is not present, the options are reset to default
118      @param options: L{SolverOptions}      :param options: `SolverOptions`
119      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
120      """      """
121      return self.__pde_p.setSolverOptions(options)      return self.__pde_p.setSolverOptions(options)
122    
# Line 124  class DarcyFlow(object): Line 124  class DarcyFlow(object):
124          """          """
125          assigns values to model parameters          assigns values to model parameters
126    
127          @param f: volumetic sources/sinks          :param f: volumetic sources/sinks
128          @type f: scalar value on the domain (e.g. L{Data})          :type f: scalar value on the domain (e.g. `Data`)
129          @param g: flux sources/sinks          :param g: flux sources/sinks
130          @type g: vector values on the domain (e.g. L{Data})          :type g: vector values on the domain (e.g. `Data`)
131          @param location_of_fixed_pressure: mask for locations where pressure is fixed          :param location_of_fixed_pressure: mask for locations where pressure is fixed
132          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
133          @param location_of_fixed_flux:  mask for locations where flux is fixed.          :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
135          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
136                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
137                               C{v} on the main diagonal is used.                               ``v`` on the main diagonal is used.
138          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})          :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
139    
140          @note: the values of parameters which are not set by calling C{setValue} are not altered.          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
142                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal                 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
143                 is along the M{x_i} axis.                 is along the *x_i* axis.
144          """          """
145          if f !=None:          if f !=None:
146             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
# Line 177  class DarcyFlow(object): Line 177  class DarcyFlow(object):
177    
178      def setTolerance(self,rtol=1e-4):      def setTolerance(self,rtol=1e-4):
179          """          """
180          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
181    
182          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183    
184          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
185    
186          @param rtol: relative tolerance for the pressure          :param rtol: relative tolerance for the pressure
187          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
188          """          """
189          if rtol<0:          if rtol<0:
190              raise ValueError,"Relative tolerance needs to be non-negative."              raise ValueError,"Relative tolerance needs to be non-negative."
# Line 193  class DarcyFlow(object): Line 193  class DarcyFlow(object):
193          """          """
194          returns the relative tolerance          returns the relative tolerance
195    
196          @return: current relative tolerance          :return: current relative tolerance
197          @rtype: C{float}          :rtype: ``float``
198          """          """
199          return self.__rtol          return self.__rtol
200    
201      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
202          """          """
203          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204    
205          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206    
207          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
208    
209          @param atol: absolute tolerance for the pressure          :param atol: absolute tolerance for the pressure
210          @type atol: non-negative C{float}          :type atol: non-negative ``float``
211          """          """
212          if atol<0:          if atol<0:
213              raise ValueError,"Absolute tolerance needs to be non-negative."              raise ValueError,"Absolute tolerance needs to be non-negative."
# Line 216  class DarcyFlow(object): Line 216  class DarcyFlow(object):
216         """         """
217         returns the absolute tolerance         returns the absolute tolerance
218                
219         @return: current absolute tolerance         :return: current absolute tolerance
220         @rtype: C{float}         :rtype: ``float``
221         """         """
222         return self.__atol         return self.__atol
223      def getSubProblemTolerance(self):      def getSubProblemTolerance(self):
224      """      """
225      Returns a suitable subtolerance      Returns a suitable subtolerance
226      @type: C{float}      @type: ``float``
227      """      """
228      return max(util.EPSILON**(0.75),self.getTolerance()**2)      return max(util.EPSILON**(0.75),self.getTolerance()**2)
229      def setSubProblemTolerance(self):      def setSubProblemTolerance(self):
# Line 244  class DarcyFlow(object): Line 244  class DarcyFlow(object):
244    
245           The iteration is terminated if the residual norm is less then self.getTolerance().           The iteration is terminated if the residual norm is less then self.getTolerance().
246    
247           @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 ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
248           @type u0: vector value on the domain (e.g. L{Data}).           :type u0: vector value on the domain (e.g. `Data`).
249           @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 ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
250           @type p0: scalar value on the domain (e.g. L{Data}).           :type p0: scalar value on the domain (e.g. `Data`).
251           @param verbose: if set some information on iteration progress are printed           :param verbose: if set some information on iteration progress are printed
252           @type verbose: C{bool}           :type verbose: ``bool``
253           @return: flux and pressure           :return: flux and pressure
254           @rtype: C{tuple} of L{Data}.           :rtype: ``tuple`` of `Data`.
255    
256           @note: The problem is solved as a least squares form           :note: The problem is solved as a least squares form
257    
258           M{(I+D^*D)u+Qp=D^*f+g}           *(I+D^*D)u+Qp=D^*f+g*
259           M{Q^*u+Q^*Qp=Q^*g}           *Q^*u+Q^*Qp=Q^*g*
260    
261           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
262           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
263    
264           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
265    
266           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
267    
268           M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure           *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
269    
270           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 *Q^*Q*). In each iteration step
271           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
272           """           """
273           self.verbose=verbose           self.verbose=verbose
274           rtol=self.getTolerance()           rtol=self.getTolerance()
275           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
276       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
277           num_corrections=0           num_corrections=0
278           converged=False           converged=False
279           p=p0           p=p0
# Line 336  class DarcyFlow(object): Line 335  class DarcyFlow(object):
335    
336      def getFlux(self,p=None, fixed_flux=Data()):      def getFlux(self,p=None, fixed_flux=Data()):
337          """          """
338          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 ``p`` where the flux is equal to ``fixed_flux``
339          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
340          Note that C{g} and C{f} are used, see L{setValue}.          Note that ``g`` and ``f`` are used, see `setValue`.
341    
342          @param p: pressure.          :param p: pressure.
343          @type p: scalar value on the domain (e.g. L{Data}).          :type p: scalar value on the domain (e.g. `Data`).
344          @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 ``location_of_fixed_flux``.
345          @type fixed_flux: vector values on the domain (e.g. L{Data}).          :type fixed_flux: vector values on the domain (e.g. `Data`).
346          @param tol: relative tolerance to be used.          :return: flux
347          @type tol: positive C{float}.          :rtype: `Data`
348          @return: flux          :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
349          @rtype: L{Data}                 for the permeability *k_{ij}*
         @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}}  
                for the permeability M{k_{ij}}  
350          """          """
351      self.setSubProblemTolerance()      self.setSubProblemTolerance()
352          g=self.__g          g=self.__g
# Line 380  class StokesProblemCartesian(Homogeneous Line 377  class StokesProblemCartesian(Homogeneous
377              sp.initialize(...)              sp.initialize(...)
378              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
379       """       """
380       def __init__(self,domain,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
381           """           """
382           initialize the Stokes Problem           initialize the Stokes Problem
383    
384           @param domain: domain of the problem. The approximation order needs to be two.           :param domain: domain of the problem. The approximation order needs to be two.
385           @type domain: L{Domain}           :type domain: `Domain`
386       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
      @type adaptSubTolerance: C{bool}  
          @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.  
387           """           """
388           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
389           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
390           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())
391           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
392            
# Line 409  class StokesProblemCartesian(Homogeneous Line 403  class StokesProblemCartesian(Homogeneous
403           """           """
404       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
405            
406       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
407       """       """
408       return self.__pde_u.getSolverOptions()       return self.__pde_u.getSolverOptions()
409       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
410           """           """
411       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
412            
413       @param options: new solver  options       :param options: new solver  options
414       @type options: L{SolverOptions}       :type options: `SolverOptions`
415       """       """
416           self.__pde_u.setSolverOptions(options)           self.__pde_u.setSolverOptions(options)
417       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
418           """           """
419       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
420       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
421       """       """
422       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
423       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
424           """           """
425       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
426       @param options: new solver  options       :param options: new solver  options
427       @type options: L{SolverOptions}       :type options: `SolverOptions`
428       """       """
429       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
430    
# Line 439  class StokesProblemCartesian(Homogeneous Line 433  class StokesProblemCartesian(Homogeneous
433       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
434       the velocity onto the function space of presure.       the velocity onto the function space of presure.
435            
436       @param options: new solver options       :param options: new solver options
437       @type options: L{SolverOptions}       :type options: `SolverOptions`
438       """       """
439       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
440       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
441           """           """
442       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
443       the velocity onto the function space of presure.       the velocity onto the function space of presure.
444            
445       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
446       """       """
447       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
      def setSubProblemTolerance(self):  
          """  
      Updates the tolerance for subproblems  
          """  
      if self.adaptSubTolerance():  
              sub_tol=self.getSubProblemTolerance()  
          self.getSolverOptionsDiv().setTolerance(sub_tol)  
          self.getSolverOptionsDiv().setAbsoluteTolerance(0.)  
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          self.getSolverOptionsVelocity().setTolerance(sub_tol)  
          self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)  
           
448    
449       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
450          """          """
451          assigns values to the model parameters          assigns values to the model parameters
452    
453          @param f: external force          :param f: external force
454          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
455          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
456          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
457          @param eta: viscosity          :param eta: viscosity
458          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
459          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
460          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
461          @param stress: initial stress          :param stress: initial stress
462      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
463          @note: All values needs to be set.          :note: All values needs to be set.
464          """          """
465          self.eta=eta          self.eta=eta
466          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
# Line 488  class StokesProblemCartesian(Homogeneous Line 469  class StokesProblemCartesian(Homogeneous
469          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
470              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
471              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
472            n=self.domain.getNormal()
473      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
474          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
475          self.__f=f          self.__f=f
476          self.__surface_stress=surface_stress          self.__surface_stress=surface_stress
477          self.__stress=stress          self.__stress=stress
478    
479       def Bv(self,v):       def Bv(self,v,tol):
480           """           """
481           returns inner product of element p and div(v)           returns inner product of element p and div(v)
482    
483           @param p: a pressure increment           :param v: a residual
484           @param v: a residual           :return: inner product of element p and div(v)
485           @return: inner product of element p and div(v)           :rtype: ``float``
486           @rtype: C{float}           """
487           """           self.__pde_proj.setValue(Y=-util.div(v)) # -???
488           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
489           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
490             out=self.__pde_proj.getSolution()
491             return out
492    
493       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
494           """           """
495           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
496    
497           @param p: a pressure increment           :param p: a pressure increment
498           @param v: a residual           :param Bv: a residual
499           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
500           @rtype: C{float}           :rtype: ``float``
501           """           """
502           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
503    
# Line 521  class StokesProblemCartesian(Homogeneous Line 505  class StokesProblemCartesian(Homogeneous
505           """           """
506           Returns inner product of p0 and p1           Returns inner product of p0 and p1
507    
508           @param p0: a pressure           :param p0: a pressure
509           @param p1: a pressure           :param p1: a pressure
510           @return: inner product of p0 and p1           :return: inner product of p0 and p1
511           @rtype: C{float}           :rtype: ``float``
512           """           """
513           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
514           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
515           return util.integrate(s0*s1)           return util.integrate(s0*s1)
516    
517       def norm_v(self,v):       def norm_v(self,v):
518           """           """
519           returns the norm of v           returns the norm of v
520    
521           @param v: a velovity           :param v: a velovity
522           @return: norm of v           :return: norm of v
523           @rtype: non-negative C{float}           :rtype: non-negative ``float``
524           """           """
525           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
526    
527       def getV(self, p, v0):       def getDV(self, p, v, tol):
528           """           """
529           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
530    
531           @param p: a pressure           :param p: a pressure
532           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
533           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
534           """           """
535           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
536         self.getSolverOptionsVelocity().setTolerance(tol)
537         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
538           if self.__stress.isEmpty():           if self.__stress.isEmpty():
539              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
540           else:           else:
541              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
542           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
543           return  out           return  out
544    
# Line 560  class StokesProblemCartesian(Homogeneous Line 546  class StokesProblemCartesian(Homogeneous
546          """          """
547          Returns Bv (overwrite).          Returns Bv (overwrite).
548    
549          @rtype: equal to the type of p          :rtype: equal to the type of p
550          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
551          """          """
552          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
553    
554       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
555           """           """
556           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
557    
558           @param p: a pressure increment           :param p: a pressure increment
559           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
560           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
561           """           """
562           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
563           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
564           return  out           return  out
565    
566       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
567           """           """
568           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
569           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
570    
571           @param v: velocity increment           :param Bv: velocity increment
572           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
573           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
574           """           """
575           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
576           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
577         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
578             out=self.__pde_prec.getSolution()
579             return out

Legend:
Removed from v.2549  
changed lines
  Added in v.2719

  ViewVC Help
Powered by ViewVC 1.1.26