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

revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC
# Line 38  from pdetools import HomogeneousSaddlePo Line 38  from pdetools import HomogeneousSaddlePo
38
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      Represents and solves the problem      solves the problem
42
43      M{u_i+k_{ij}*p_{,j} = g_i}      M{u_i+k_{ij}*p_{,j} = g_i}

44      M{u_{i,i} = f}      M{u_{i,i} = f}
45
46      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
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):      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
55          @type domain: L{Domain}          @type domain: L{Domain}
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          """          """
71          Assigns values to model parameters.          assigns values to model parameters
72
73          @param f: volumetric sources/sinks          @param f: volumetic sources/sinks
74          @type f: scalar value on the domain, e.g. L{Data}          @type f: scalar value on the domain (e.g. L{Data})
75          @param g: flux sources/sinks          @param g: flux sources/sinks
76          @type g: vector value on the domain, e.g. L{Data}          @type g: vector values on the domain (e.g. L{Data})
77          @param location_of_fixed_pressure: mask for locations where pressure is fixed          @param location_of_fixed_pressure: mask for locations where pressure is fixed
78          @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. L{Data})
79          @param location_of_fixed_flux: mask for locations where flux is fixed          @param location_of_fixed_flux:  mask for locations where flux is fixed.
80          @type location_of_fixed_flux: vector value on the domain (e.g. L{Data})          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
81          @param permeability: permeability tensor. If scalar C{s} is given the          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
82                               tensor with C{s} on the main diagonal is used. If                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
83                               vector C{v} is given the tensor with C{v} on the                               C{v} on the main diagonal is used.
84                               main diagonal is used.          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
85          @type permeability: scalar, vector or tensor values on the domain, e.g.
86                              L{Data}          @note: the values of parameters which are not set by calling C{setValue} are not altered.
87            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
88          @note: the values of parameters which are not set by calling                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
89                 C{setValue} are not altered                 is along the M{x_i} axis.
@note: at any point on the boundary of the domain the pressure
(C{location_of_fixed_pressure}) >0 or the normal component of
the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
the normal is along the M{x_i} axis.
90          """          """
91          if f !=None:          if f !=None:
92             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
93             if f.isEmpty():             if f.isEmpty():
94                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
95             else:             else:
96                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
97             self.f=f             self.f=f
98          if g !=None:          if g !=None:
99             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
100             if g.isEmpty():             if g.isEmpty():
101               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 126  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}.          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}).
129          The flux is equal to C{fixed_flux} on locations where          Note that C{g} and C{f} are used, see L{setValue}.
130          C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}
131          and C{f} are used.          @param p: pressure.
132            @type p: scalar value on the domain (e.g. L{Data}).
133          @param p: pressure          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134          @type p: scalar value on the domain, e.g. L{Data}          @type fixed_flux: vector values on the domain (e.g. L{Data}).
135          @param fixed_flux: flux on the locations of the domain marked by          @param tol: relative tolerance to be used.
136                             C{location_of_fixed_flux}          @type tol: positive C{float}.
137          @type fixed_flux: vector values on the domain, e.g. L{Data}          @return: flux
@param tol: relative tolerance to be used
@type tol: positive float
@return: flux
138          @rtype: L{Data}          @rtype: L{Data}
139          @note: the method uses the least squares solution          @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                 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator                 for the permeability M{k_{ij}}
and M{(Qp)_i=k_{ij}p_{,j}} 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          Solves the problem.          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
The iteration is terminated if the error in the pressure is less than
M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
side (see escript user's guide for details).

@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.
@type u0: vector value on the domain, e.g. L{Data}
@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.
@type p0: scalar value on the domain, e.g. L{Data}
196          @param atol: absolute tolerance for the pressure          @param atol: absolute tolerance for the pressure
197          @type atol: non-negative C{float}          @type atol: non-negative C{float}
198          @param rtol: relative tolerance for the pressure          @param rtol: relative tolerance for the pressure
199          @type rtol: non-negative C{float}          @type rtol: non-negative C{float}
200          @param sub_rtol: tolerance to be used in the sub iteration. It is          @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                           recommended that M{sub_rtol<rtol*5.e-3}          L{getPressure} method or use  the value from a previous time step.
202          @type sub_rtol: positive-negative C{float}          @type p_ref: scalar value on the domain (e.g. L{Data}).
203          @param verbose: if True information on iteration progress is printed          @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          @type verbose: C{bool}          L{getFlux} method or use  the value from a previous time step.
205          @param show_details: if True information on the sub-iteration process          @type v_ref: vector-valued on the domain (e.g. L{Data}).
206                               is printed          @return: used absolute tolerance.
207          @type show_details: C{bool}          @rtype: positive C{float}
208          @return: flux and pressure          """
209          @rtype: C{tuple} of L{Data}          g=self.__g
210            if not v_ref == None:
211          @note: the problem is solved in a least squares formulation:             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212            else:
213          M{(I+D^*D)u+Qp=D^*f+g}             f1=util.integrate(util.length(util.interpolate(g))**2)
214            if not p_ref == None:
215          M{Q^*u+Q^*Qp=Q^*g}             f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
216            else:
217          where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the             f2=0
218          permeability M{k_{ij}}. We eliminate the flux from the problem by          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219          setting          if self.__ATOL<=0:
220               raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221          M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}          return self.__ATOL
222
223          from the first equation. Inserted into the second equation we get      def getTolerance(self):
224            """
225          M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}          returns the current tolerance.
226          on C{location_of_fixed_pressure}
227            @return: used absolute tolerance.
228          which is solved using the PCG method (precondition is M{Q^*Q}).          @rtype: positive C{float}
229          In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}          """
230          need to be solved using a sub-iteration scheme.          if self.__ATOL==None:
231          """             raise ValueError,"no tolerance is defined."
232          self.verbose=verbose          return self.__ATOL
233          self.show_details= show_details and self.verbose
234          self.__pde_v.setTolerance(sub_rtol)      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235          self.__pde_p.setTolerance(sub_rtol)           """
236          u2=u0*self.__pde_v.getCoefficient("q")           solves the problem.
237          #
238          # first the reference velocity is calculated from           The iteration is terminated if the residual norm is less then self.getTolerance().
239          #
240          #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)           @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}).
242          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)           @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          u_ref=self.__pde_v.getSolution(verbose=show_details)           @type p0: scalar value on the domain (e.g. L{Data}).
244          if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245          self.__pde_v.setValue(r=Data())           @type sub_rtol: positive-negative C{float}
246          #           @param verbose: if set some information on iteration progress are printed
247          #   and then we calculate a reference pressure           @type verbose: C{bool}
248          #           @param show_details:  if set information on the subiteration process are printed.
249          #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)           @type show_details: C{bool}
250          #           @return: flux and pressure
251          self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)           @rtype: C{tuple} of L{Data}.
252          p_ref=self.__pde_p.getSolution(verbose=self.show_details)
253          if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)           @note: The problem is solved as a least squares form
254          self.__pde_p.setValue(r=Data())
255          #           M{(I+D^*D)u+Qp=D^*f+g}
256          #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref           M{Q^*u+Q^*Qp=Q^*g}
257          #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
258          #           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
259          #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du           We eliminate the flux form the problem by setting
260          #
261          #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
262          #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
263          #           form the first equation. Inserted into the second equation we get
264          #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
265          #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)           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
266          #
267          #  we use the (du,Qdp) to represent the resudual           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
268          #  Q^*Q is a preconditioner           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
269          #           """
270          #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>           self.verbose=verbose
271          #           self.show_details= show_details and self.verbose
272          Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))           self.__pde_v.setTolerance(sub_rtol)
273          norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))           self.__pde_p.setTolerance(sub_rtol)
274          ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)           ATOL=self.getTolerance()
275          if not ATOL>0:           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276              raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)           #########################################################################################################################
277          if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)           #
278          #           #   we solve:
279          #   caclulate the initial residual           #
280          #           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
281          self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())           #
282          du=self.__pde_v.getSolution(verbose=show_details)           #   residual is
283          r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)           #
284          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)           #    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)
285          util.saveVTK("d.vtu",p=dp,p_ref=p_ref)           #
286          return u_ref+r[1],dp           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
287             #
288      def __Aprod_PCG(self,p):           #    we use (g - Qp, v) to represent the residual. not that
289          if self.show_details: print "DarcyFlux: Applying operator"           #
290          Qp=util.tensor_mult(self.__permeability,util.grad(p))           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291          self.__pde_v.setValue(Y=Qp,X=Data())           #
292          w=self.__pde_v.getSolution(verbose=self.show_details)           #   while the initial residual is
293          return ArithmeticTuple(-Qp,w)           #
294             #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295             #
296             d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297             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             if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300             self.__pde_v.setValue(r=Data())
301             # start CG
302             r=ArithmeticTuple(d0, v00)
303             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             return r[1],p
305
306        def __Aprod_PCG(self,dp):
307              if self.show_details: print "DarcyFlux: Applying operator"
308              #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
309              mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
310              self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311              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        Represents and solves the problem        solves
329
330        M{-(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
332
333        M{u_{i,i}=0} and M{u=0} where C{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
335
336        M{eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j}        if surface_stress is not given 0 is assumed.
337
338        If C{surface_stress} is not given 0 is assumed.        typical usage:
339
340        Typical usage::              sp=StokesProblemCartesian(domain)
341                sp.setTolerance()
342            sp = StokesProblemCartesian(domain)              sp.initialize(...)
343            sp.setTolerance()              v,p=sp.solve(v0,p0)
sp.initialize(...)
v,p = sp.solve(v0,p0)
344        """        """
345        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
346           """           """
347           Initializes the Stokes Problem.           initialize the Stokes Problem
348
349           @param domain: domain of the problem. The approximation order needs           @param domain: domain of the problem. The approximation order needs to be two.
to be two.
350           @type domain: L{Domain}           @type domain: L{Domain}
351           @warning: The approximation order needs to be two otherwise you may           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
see oscillations in the pressure.
352           """           """
353           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
354           self.domain=domain           self.domain=domain
# Line 312  class StokesProblemCartesian(Homogeneous Line 357  class StokesProblemCartesian(Homogeneous
357           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
358           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
359           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
360
361           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
362           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
363           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
# Line 324  class StokesProblemCartesian(Homogeneous Line 369  class StokesProblemCartesian(Homogeneous
369           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
370
371        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()):
372           """          """
373           Assigns values to the model parameters.          assigns values to the model parameters
374
375           @param f: external force          @param f: external force
376           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
377           @param fixed_u_mask: mask of locations with fixed velocity          @param fixed_u_mask: mask of locations with fixed velocity.
378           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
379                               or similar          @param eta: viscosity
380           @param eta: viscosity          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
381           @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar          @param surface_stress: normal surface stress
382           @param surface_stress: normal surface stress          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
383           @type surface_stress: L{Vector} object on L{FunctionSpace},          @param stress: initial stress
384                                 L{FunctionOnBoundary} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
385           @param stress: initial stress          @note: All values needs to be set.
386           @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or
387                         similar          """
388           @note: All values need to be set.          self.eta=eta
389           """          A =self.__pde_u.createCoefficient("A")
390           self.eta=eta      self.__pde_u.setValue(A=Data())
391           A =self.__pde_u.createCoefficient("A")          for i in range(self.domain.getDim()):
392           self.__pde_u.setValue(A=Data())          for j in range(self.domain.getDim()):
393           for i in range(self.domain.getDim()):              A[i,j,j,i] += 1.
394               for j in range(self.domain.getDim()):              A[i,j,i,j] += 1.
395                   A[i,j,j,i] += 1.      self.__pde_prec.setValue(D=1/self.eta)
396                   A[i,j,i,j] += 1.          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
397           self.__pde_prec.setValue(D=1/self.eta)          self.__stress=stress
self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
self.__stress=stress
398
399        def B(self,v):        def B(self,v):
400           """          """
401           Returns M{div(v)}.          returns div(v)
402           @return: M{div(v)}          @rtype: equal to the type of p
@rtype: equal to the type of p
403
404           @note: Boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
405           """          """
406           if self.show_details: print "apply divergence:"          if self.show_details: print "apply divergence:"
407           self.__pde_proj.setValue(Y=-util.div(v))          self.__pde_proj.setValue(Y=-util.div(v))
408           self.__pde_proj.setTolerance(self.getSubProblemTolerance())          self.__pde_proj.setTolerance(self.getSubProblemTolerance())
409           return self.__pde_proj.getSolution(verbose=self.show_details)          return self.__pde_proj.getSolution(verbose=self.show_details)
410
411        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
412           """           """
413           Returns inner product of element p and Bv (overwrite).           returns inner product of element p and Bv  (overwrite)
414
415           @type p: equal to the type of p           @type p: equal to the type of p
416           @type Bv: equal to the type of result of operator B           @type Bv: equal to the type of result of operator B
417           @return: inner product of p and Bv           @rtype: C{float}
418
419           @rtype: equal to the type of p           @rtype: equal to the type of p
420           """           """
421           s0=util.interpolate(p,Function(self.domain))           s0=util.interpolate(p,Function(self.domain))
# Line 381  class StokesProblemCartesian(Homogeneous Line 424  class StokesProblemCartesian(Homogeneous
424
425        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
426           """           """
427           Returns inner product of element p0 and p1 (overwrite).           returns inner product of element p0 and p1  (overwrite)
428
429           @type p0: equal to the type of p           @type p0: equal to the type of p
430           @type p1: equal to the type of p           @type p1: equal to the type of p
431           @return: inner product of p0 and p1           @rtype: C{float}
432
433           @rtype: equal to the type of p           @rtype: equal to the type of p
434           """           """
435           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
# Line 394  class StokesProblemCartesian(Homogeneous Line 438  class StokesProblemCartesian(Homogeneous
438
439        def inner_v(self,v0,v1):        def inner_v(self,v0,v1):
440           """           """
441           Returns inner product of two elements v0 and v1 (overwrite).           returns inner product of two element v0 and v1  (overwrite)
442
443           @type v0: equal to the type of v           @type v0: equal to the type of v
444           @type v1: equal to the type of v           @type v1: equal to the type of v
445           @return: inner product of v0 and v1           @rtype: C{float}
446
447           @rtype: equal to the type of v           @rtype: equal to the type of v
448           """           """
449           gv0=util.grad(v0)       gv0=util.grad(v0)
450           gv1=util.grad(v1)       gv1=util.grad(v1)
451           return util.integrate(util.inner(gv0,gv1))           return util.integrate(util.inner(gv0,gv1))
452
453        def solve_A(self,u,p):        def solve_A(self,u,p):
454           """           """
455           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
456           """           """
457           if self.show_details: print "solve for velocity:"           if self.show_details: print "solve for velocity:"
458           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
# Line 416  class StokesProblemCartesian(Homogeneous Line 461  class StokesProblemCartesian(Homogeneous
461           else:           else:
462              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-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
463           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
464           return out           return  out
465
466        def solve_prec(self,p):        def solve_prec(self,p):
"""
Applies the preconditioner.
"""
467           if self.show_details: print "apply preconditioner:"           if self.show_details: print "apply preconditioner:"
468           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
469           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
470           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
471           return q           return q

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

 ViewVC Help Powered by ViewVC 1.1.26