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

revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC revision 2415 by gross, Wed May 13 02:48:39 2009 UTC
# Line 16  http://www.uq.edu.au/esscc Line 16  http://www.uq.edu.au/esscc
16  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
17  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
18  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
19  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
20
21  """  """
22  Some models for flow  Some models for flow
# Line 34  __author__="Lutz Gross, l.gross@uq.edu.a Line 34  __author__="Lutz Gross, l.gross@uq.edu.a
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      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 the permeability,      where M{p} represents the pressure and M{u} the Darcy flux. M{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):      def __init__(self, domain, weight=None, 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            if weight == None:
59               s=self.domain.getSize()
60               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61            else:
62               self.__l=weight
63          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
64          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))          if useReduced: self.__pde_v.setReducedOrderOn()
65          self.__pde_v.setSymmetryOn()          self.__pde_v.setSymmetryOn()
66            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68          self.__pde_p=LinearSinglePDE(domain)          self.__pde_p=LinearSinglePDE(domain)
69          self.__pde_p.setSymmetryOn()          self.__pde_p.setSymmetryOn()
70            if useReduced: self.__pde_p.setReducedOrderOn()
71          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73            self.setTolerance()
74            self.setAbsoluteTolerance()
75            self.setSubProblemTolerance()
76
77      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):
78          """          """
# Line 75  class DarcyFlow(object): Line 86  class DarcyFlow(object):
86          @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})
87          @param location_of_fixed_flux:  mask for locations where flux is fixed.          @param location_of_fixed_flux:  mask for locations where flux is fixed.
88          @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. L{Data})
89          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91                               C{v} on the main diagonal is used.                               C{v} on the main diagonal is used.
92          @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. L{Data})
93
# Line 85  class DarcyFlow(object): Line 96  class DarcyFlow(object):
96                 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 (C{location_of_fixed_flux[i]>0} if direction of the normal
97                 is along the M{x_i} axis.                 is along the M{x_i} axis.
98          """          """
99          if f !=None:          if f !=None:
100             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101             if f.isEmpty():             if f.isEmpty():
102                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103             else:             else:
104                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
105             self.f=f             self.__f=f
106          if g !=None:            if g !=None:
107             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108             if g.isEmpty():             if g.isEmpty():
109               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 118  class DarcyFlow(object): Line 129  class DarcyFlow(object):
129             self.__permeability=perm             self.__permeability=perm
130             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))
131
132        def setTolerance(self,rtol=1e-4):
133            """
134            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135
136            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137
138            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}}.
139
140      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8):          @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142          """          """
143          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          if rtol<0:
144          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).              raise ValueError,"Relative tolerance needs to be non-negative."
145          Note that C{g} and C{f} are used, L{setValue}.          self.__rtol=rtol
146                def getTolerance(self):
147          @param p: pressure.          """
148          @type p: scalar value on the domain (e.g. L{Data}).          returns the relative tolerance
149          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
150          @type fixed_flux: vector values on the domain (e.g. L{Data}).          @return: current relative tolerance
151          @param tol: relative tolerance to be used.          @rtype: C{float}
152          @type tol: positive float.          """
153          @return: flux          return self.__rtol
154          @rtype: L{Data}
155          @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}}      def setAbsoluteTolerance(self,atol=0.):
for the permeability M{k_{ij}}
156          """          """
157          self.__pde_v.setTolerance(tol)          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=boundary_flux)
159          return self.__pde_v.getSolution()          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160
161            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}}.
162
163            @param atol: absolute tolerance for the pressure
164            @type atol: non-negative C{float}
165            """
166            if atol<0:
167                raise ValueError,"Absolute tolerance needs to be non-negative."
168            self.__atol=atol
169        def getAbsoluteTolerance(self):
170           """
171           returns the absolute tolerance
172
173           @return: current absolute tolerance
174           @rtype: C{float}
175           """
176           return self.__atol
177
178        def setSubProblemTolerance(self,rtol=None):
179             """
180             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181             C{self.getTolerance()**2} is used.
182
183             @param rtol: relative tolerence
184             @type rtol: positive C{float}
185             """
186             if rtol == None:
187                  if self.getTolerance()<=0.:
188                      raise ValueError,"A positive relative tolerance must be set."
189                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190             else:
191                 if rtol<=0:
192                     raise ValueError,"sub-problem tolerance must be positive."
193                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194
195      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def getSubProblemTolerance(self):
196           """           """
197             Returns the subproblem reduction factor.
198
199             @return: subproblem reduction factor
200             @rtype: C{float}
201             """
202             return self.__sub_tol
203
204        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205             """
206           solves the problem.           solves the problem.
207
208           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).
209
210           @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.
211           @type u0: vector value on the domain (e.g. L{Data}).           @type u0: vector value on the domain (e.g. L{Data}).
212           @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.
213           @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}
@param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
@type sub_rtol: positive-negative C{float}
214           @param verbose: if set some information on iteration progress are printed           @param verbose: if set some information on iteration progress are printed
215           @type verbose: C{bool}           @type verbose: C{bool}
216           @param show_details:  if set information on the subiteration process are printed.           @param show_details:  if set information on the subiteration process are printed.
217           @type show_details: C{bool}           @type show_details: C{bool}
218           @return: flux and pressure           @return: flux and pressure
219           @rtype: C{tuple} of L{Data}.           @rtype: C{tuple} of L{Data}.
220
221           @note: The problem is solved as a least squares form           @note: The problem is solved as a least squares form
222
223           M{(I+D^*D)u+Qp=D^*f+g}           M{(I+D^*D)u+Qp=D^*f+g}
224           M{Q^*u+Q^*Qp=Q^*g}           M{Q^*u+Q^*Qp=Q^*g}
225
226           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227           We eliminate the flux form the problem by setting           We eliminate the flux form the problem by setting
228
229           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
230
231           form the first equation. Inserted into the second equation we get           form the first equation. Inserted into the second equation we get
232
233           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           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
234
235           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 M{Q^*Q}). In each iteration step
236           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237           """           """
238           self.verbose=verbose           self.verbose=verbose
239           self.show_details= show_details and self.verbose           self.show_details= show_details and self.verbose
240           self.__pde_v.setTolerance(sub_rtol)           rtol=self.getTolerance()
241           self.__pde_p.setTolerance(sub_rtol)           atol=self.getAbsoluteTolerance()
242           p2=p0*self.__pde_p.getCoefficient("q")           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243           u2=u0*self.__pde_v.getCoefficient("q")
244           g=self.__g-u2-util.tensor_mult(self.__permeability,util.grad(p2))           num_corrections=0
245           f=self.__f-util.div(u2)           converged=False
246           self.__pde_v.setValue(Y=g, X=f*util.kronecker(self.domain), r=Data())           p=p0
247           dv=self.__pde_v.getSolution(verbose=show_details)           norm_r=None
248           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-dv))           while not converged:
249           self.__pde_p.setValue(r=Data())                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250           dp=self.__pde_p.getSolution(verbose=self.show_details)                 Qp=self.__Q(p)
251           norm_rhs=self.__inner_PCG(dp,ArithmeticTuple(g,dv))                 norm_v=self.__L2(v)
252           if norm_rhs<0:                 norm_Qp=self.__L2(Qp)
253               raise NegativeNorm,"negative norm. Maybe the sub-tolerance is too large."                 if norm_v == 0.:
254           ATOL=util.sqrt(norm_rhs)*rtol +atol                    if norm_Qp == 0.:
255           if not ATOL>0:                       return v,p
256               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol)                    else:
257           rhs=ArithmeticTuple(g,dv)                      fac=norm_Qp
258           dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0, verbose=self.verbose, initial_guess=False)                 else:
259           return u2+r[1],p2+dp                    if norm_Qp == 0.:
260                                fac=norm_v
261      def __Aprod_PCG(self,p):                    else:
262            if self.show_details: print "DarcyFlux: Applying operator"                      fac=2./(1./norm_v+1./norm_Qp)
263            Qp=util.tensor_mult(self.__permeability,util.grad(p))                 ATOL=(atol+rtol*fac)
264            self.__pde_v.setValue(Y=Qp,X=Data())                 if self.verbose:
265            w=self.__pde_v.getSolution(verbose=self.show_details)                      print "DarcyFlux: L2 norm of v = %e."%norm_v
266            return ArithmeticTuple(Qp,w)                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268                   if norm_r == None or norm_r>ATOL:
269                       if num_corrections>max_num_corrections:
270                             raise ValueError,"maximum number of correction steps reached."
271                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272                       num_corrections+=1
273                   else:
274                       converged=True
275             return v,p
276        def __L2(self,v):
277             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278
279        def __Q(self,p):
280              return util.tensor_mult(self.__permeability,util.grad(p))
281
282        def __Aprod(self,dp):
283              self.__pde_v.setTolerance(self.getSubProblemTolerance())
284              if self.show_details: print "DarcyFlux: Applying operator"
285              Qdp=self.__Q(dp)
286              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288              # self.__pde_v.getOperator().saveMM("proj.mm")
289              return Qdp+du
290        def __inner_GMRES(self,r,s):
291             return util.integrate(util.inner(r,s))
292
293      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
294           a=util.tensor_mult(self.__permeability,util.grad(p))           return util.integrate(util.inner(self.__Q(p), r))
return util.integrate(util.inner(a,r[0]-r[1]))
295
296      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
297              self.__pde_p.setTolerance(self.getSubProblemTolerance())
298            if self.show_details: print "DarcyFlux: Applying preconditioner"            if self.show_details: print "DarcyFlux: Applying preconditioner"
299            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), Y=Data(), r=Data())
300            return self.__pde_p.getSolution(verbose=self.show_details)            # self.__pde_p.getOperator().saveMM("prec.mm")
301              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302
303        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304            """
305            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
306            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
307            Note that C{g} and C{f} are used, see L{setValue}.
308
309            @param p: pressure.
310            @type p: scalar value on the domain (e.g. L{Data}).
311            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
312            @type fixed_flux: vector values on the domain (e.g. L{Data}).
313            @param tol: relative tolerance to be used.
314            @type tol: positive C{float}.
315            @return: flux
316            @rtype: L{Data}
317            @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}}
318                   for the permeability M{k_{ij}}
319            """
320            self.__pde_v.setTolerance(self.getSubProblemTolerance())
321            g=self.__g
322            f=self.__f
323            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
324            if p == None:
325               self.__pde_v.setValue(Y=g)
326            else:
327               self.__pde_v.setValue(Y=g-self.__Q(p))
328            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329
330  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
331        """       """
332        solves       solves
333
334            -(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}
335                  u_{i,i}=0                  u_{i,i}=0
# Line 230  class StokesProblemCartesian(Homogeneous Line 337  class StokesProblemCartesian(Homogeneous
337            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
338            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
339
340        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
341
342        typical usage:       typical usage:
343
344              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
345              sp.setTolerance()              sp.setTolerance()
346              sp.initialize(...)              sp.initialize(...)
347              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
348        """       """
349        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
350           """           """
351           initialize the Stokes Problem           initialize the Stokes Problem
352
# Line 253  class StokesProblemCartesian(Homogeneous Line 360  class StokesProblemCartesian(Homogeneous
360           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())
361           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
362           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364
365           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
366           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
367           self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
369
370           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
371           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
372         self.__pde_proj.setValue(D=1)
373           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
self.__pde_proj.setValue(D=1.)
374
375        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
376         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
377          """          """
378          assigns values to the model parameters          assigns values to the model parameters
379
# Line 280  class StokesProblemCartesian(Homogeneous Line 388  class StokesProblemCartesian(Homogeneous
388          @param stress: initial stress          @param stress: initial stress
389      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
390          @note: All values needs to be set.          @note: All values needs to be set.
391
392          """          """
393          self.eta=eta          self.eta=eta
394          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
395      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
396          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
397          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
398              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
399              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
400      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
401          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)
402            self.__f=f
403            self.__surface_stress=surface_stress
404          self.__stress=stress          self.__stress=stress
405
406        def B(self,v):       def inner_pBv(self,p,v):
407          """           """
408          returns div(v)           returns inner product of element p and div(v)
@rtype: equal to the type of p
409
410          @note: boundary conditions on p should be zero!           @param p: a pressure increment
411          """           @param v: a residual
412          if self.show_details: print "apply divergence:"           @return: inner product of element p and div(v)
self.__pde_proj.setValue(Y=-util.div(v))
self.__pde_proj.setTolerance(self.getSubProblemTolerance())
return self.__pde_proj.getSolution(verbose=self.show_details)

def inner_pBv(self,p,Bv):
"""
returns inner product of element p and Bv  (overwrite)

@type p: equal to the type of p
@type Bv: equal to the type of result of operator B
413           @rtype: C{float}           @rtype: C{float}

@rtype: equal to the type of p
414           """           """
415           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)
416
417        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
418           """           """
419           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}
420
421           @rtype: equal to the type of p           @param p0: a pressure
422             @param p1: a pressure
423             @return: inner product of p0 and p1
424             @rtype: C{float}
425           """           """
426           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
427           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1/self.eta,Function(self.domain))
428           return util.integrate(s0*s1)           return util.integrate(s0*s1)
429
430        def inner_v(self,v0,v1):       def norm_v(self,v):
431           """           """
432           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}
433
434           @rtype: equal to the type of v           @param v: a velovity
435             @return: norm of v
436             @rtype: non-negative C{float}
437           """           """
438       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))
439
440        def solve_A(self,u,p):       def getV(self, p, v0):
441           """           """
442           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
443
444             @param p: a pressure
445             @param v0: a initial guess for the value v to return.
446             @return: v given as M{v= A^{-1} (f-B^*p)}
447           """           """
if self.show_details: print "solve for velocity:"
448           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
449             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
450           if self.__stress.isEmpty():           if self.__stress.isEmpty():
451              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))
452           else:           else:
453              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))
454             out=self.__pde_u.getSolution(verbose=self.show_details)
455             return  out
456
457
458             raise NotImplementedError,"no v calculation implemented."
459
460
461         def norm_Bv(self,v):
462            """
463            Returns Bv (overwrite).
464
465            @rtype: equal to the type of p
466            @note: boundary conditions on p should be zero!
467            """
468            self.__pde_proj.setValue(Y=util.div(v))
469            self.__pde_prec.setTolerance(self.getSubProblemTolerance())
470            return util.sqrt(util.integrate(util.interpolate(self.__pde_proj.getSolution(),Function(self.domain))**2))
471
472         def solve_AinvBt(self,p):
473             """
474             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
475
476             @param p: a pressure increment
477             @return: the solution of M{Av=B^*p}
478             @note: boundary conditions on v should be zero!
479             """
480             self.__pde_u.setTolerance(self.getSubProblemTolerance())
481             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
482           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
483           return  out           return  out
484
485        def solve_prec(self,p):       def solve_precB(self,v):
486           if self.show_details: print "apply preconditioner:"           """
487             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
488             with accuracy L{self.getSubProblemTolerance()} (overwrite).
489
490             @param v: velocity increment
491             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
492             @note: boundary conditions on p are zero.
493             """
494             self.__pde_prec.setValue(Y=-util.div(v))
495           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
496           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.2100 changed lines Added in v.2415

 ViewVC Help Powered by ViewVC 1.1.26