# Diff of /branches/restext/escript/py_src/flows.py

revision 2615 by jfenwick, Wed Aug 19 05:44:15 2009 UTC revision 2616 by jfenwick, Wed Aug 19 06:33:08 2009 UTC
# 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 *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,      where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
47
# Line 181  class DarcyFlow(object): Line 181  class DarcyFlow(object):
181
182          *|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 ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|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 ``float``          :type rtol: non-negative ``float``
# Line 204  class DarcyFlow(object): Line 204  class DarcyFlow(object):
204
205          *|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 ``rtol`` is an absolut tolerance (see `setTolerance`), *|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 ``float``          :type atol: non-negative ``float``
# Line 258  class DarcyFlow(object): Line 258  class DarcyFlow(object):
258           *(I+D^*D)u+Qp=D^*f+g*           *(I+D^*D)u+Qp=D^*f+g*
259           *Q^*u+Q^*Qp=Q^*g*           *Q^*u+Q^*Qp=Q^*g*
260
261           where *D* is the *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 *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 *I+D^*D* and with *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.
# Line 344  class DarcyFlow(object): Line 344  class DarcyFlow(object):
344          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `Data`).
345          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
346          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type fixed_flux: vector values on the domain (e.g. `Data`).
:param tol: relative tolerance to be used.
:type tol: positive ``float``.
347          :return: flux          :return: flux
348          :rtype: `Data`          :rtype: `Data`
349          :note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where *D* is the *div* operator and M{(Qp)_i=k_{ij}p_{,j}}          :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}*
350                 for the permeability M{k_{ij}}                 for the permeability *k_{ij}*
351          """          """
352      self.setSubProblemTolerance()      self.setSubProblemTolerance()
353          g=self.__g          g=self.__g
# Line 476  class StokesProblemCartesian(Homogeneous Line 474  class StokesProblemCartesian(Homogeneous
474          :param eta: viscosity          :param eta: viscosity
475          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
476          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
477          :type eta: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
478          :param stress: initial stress          :param stress: initial stress
479      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
480          :note: All values needs to be set.          :note: All values needs to be set.
# Line 498  class StokesProblemCartesian(Homogeneous Line 496  class StokesProblemCartesian(Homogeneous
496           """           """
497           returns inner product of element p and div(v)           returns inner product of element p and div(v)
498
:param p: a pressure increment
499           :param v: a residual           :param v: a residual
500           :return: inner product of element p and div(v)           :return: inner product of element p and div(v)
501           :rtype: ``float``           :rtype: ``float``
# Line 511  class StokesProblemCartesian(Homogeneous Line 508  class StokesProblemCartesian(Homogeneous
508           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
509
510           :param p: a pressure increment           :param p: a pressure increment
511           :param v: a residual           :param Bv: a residual
512           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
513           :rtype: ``float``           :rtype: ``float``
514           """           """
# Line 546  class StokesProblemCartesian(Homogeneous Line 543  class StokesProblemCartesian(Homogeneous
543
544           :param p: a pressure           :param p: a pressure
545           :param v0: a initial guess for the value v to return.           :param v0: a initial guess for the value v to return.
546           :return: v given as M{v= A^{-1} (f-B^*p)}           :return: v given as *v= A^{-1} (f-B^*p)*
547           """           """
548           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
549           if self.__stress.isEmpty():           if self.__stress.isEmpty():
# Line 579  class StokesProblemCartesian(Homogeneous Line 576  class StokesProblemCartesian(Homogeneous
576
577       def solve_prec(self,Bv):       def solve_prec(self,Bv):
578           """           """
579           applies preconditioner for for M{BA^{-1}B^*} to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
580           with accuracy `self.getSubProblemTolerance()`           with accuracy `self.getSubProblemTolerance()`
581
582           :param v: velocity increment           :param Bv: velocity increment
583           :return: *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^ * )*
584           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
585           """           """
586           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)

Legend:
 Removed from v.2615 changed lines Added in v.2616