/[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 2168 by gross, Mon Dec 15 05:09:02 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC
# Line 38  from pdetools import HomogeneousSaddlePo Line 38  from pdetools import HomogeneousSaddlePo
38    
39  class DarcyFlow(object):  class DarcyFlow(object):
40      """      """
41      solves the problem      Represents and solves the problem
42    
43      M{u_i+k_{ij}*p_{,j} = g_i}      M{u_i+k_{ij}*p_{,j} = g_i}
44    
45      M{u_{i,i} = f}      M{u_{i,i} = f}
46    
47      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
48        the permeability.
49    
50      @note: The problem is solved in a least squares formulation.      @note: The problem is solved in a least squares formulation.
51      """      """
52    
53      def __init__(self, domain):      def __init__(self, domain):
54          """          """
55          initializes the Darcy flux problem          Initializes the Darcy flux problem.
56    
57          @param domain: domain of the problem          @param domain: domain of the problem
58          @type domain: L{Domain}          @type domain: L{Domain}
59          """          """
# Line 65  class DarcyFlow(object): Line 68  class DarcyFlow(object):
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: volumetic sources/sinks          @param f: volumetric 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 values on the domain (e.g. L{Data})          @type g: vector value 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 values on the domain (e.g. L{Data})          @type location_of_fixed_flux: vector value on the domain (e.g. L{Data})
81          @param permeability: permeability tensor. If scalar C{s} is given the tensor with          @param permeability: permeability tensor. If scalar C{s} is given the
82                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with                               tensor with C{s} on the main diagonal is used. If
83                               C{v} on the main diagonal is used.                               vector C{v} is given the tensor with C{v} on the
84          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})                               main diagonal is used.
85            @type permeability: scalar, vector or tensor values on the domain, e.g.
86          @note: the values of parameters which are not set by calling C{setValue} are not altered.                              L{Data}
87          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)  
88                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal          @note: the values of parameters which are not set by calling
89                 is along the M{x_i} axis.                 C{setValue} are not altered
90            @note: at any point on the boundary of the domain the pressure
91                   (C{location_of_fixed_pressure}) >0 or the normal component of
92                   the flux (C{location_of_fixed_flux[i]}) >0 if the direction of
93                   the normal is along the M{x_i} axis.
94          """          """
95          if f !=None:          if f !=None:
96             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
97             if f.isEmpty():             if f.isEmpty():
98                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
99             else:             else:
100                 if f.getRank()>0: raise ValueError,"illegal rank of f."                 if f.getRank()>0: raise ValueError,"illegal rank of f."
101             self.f=f             self.f=f
102          if g !=None:            if g !=None:
103             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
104             if g.isEmpty():             if g.isEmpty():
105               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
# Line 121  class DarcyFlow(object): Line 128  class DarcyFlow(object):
128    
129      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):      def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
130          """          """
131          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}          Returns the flux for a given pressure C{p}.
132          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
133          Note that C{g} and C{f} are used, L{setValue}.          The flux is equal to C{fixed_flux} on locations where
134                    C{location_of_fixed_flux} is positive (see L{setValue}). Note that C{g}
135          @param p: pressure.          and C{f} are used.
136          @type p: scalar value on the domain (e.g. L{Data}).  
137          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.          @param p: pressure
138          @type fixed_flux: vector values on the domain (e.g. L{Data}).          @type p: scalar value on the domain, e.g. L{Data}
139          @param tol: relative tolerance to be used.          @param fixed_flux: flux on the locations of the domain marked by
140          @type tol: positive float.                             C{location_of_fixed_flux}
141          @return: flux          @type fixed_flux: vector values on the domain, e.g. L{Data}
142            @param tol: relative tolerance to be used
143            @type tol: positive float
144            @return: flux
145          @rtype: L{Data}          @rtype: L{Data}
146          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}          @note: the method uses the least squares solution
147                 for the permeability M{k_{ij}}                 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator
148                   and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
149          """          """
150          self.__pde_v.setTolerance(tol)          self.__pde_v.setTolerance(tol)
151          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
152          return self.__pde_v.getSolution(verbose=show_details)          return self.__pde_v.getSolution(verbose=show_details)
153    
154      def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):      def solve(self, u0, p0, atol=0, rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
155           """          """
156           solves the problem.          Solves the problem.
157    
158           The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where          The iteration is terminated if the error in the pressure is less than
159           C{|q|} denotes the norm of the right hand side (see escript user's guide for details).          M{rtol * |q| + atol} where M{|q|} denotes the norm of the right hand
160            side (see escript user's guide for details).
161           @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.  
162           @type u0: vector value on the domain (e.g. L{Data}).          @param u0: initial guess for the flux. At locations in the domain
163           @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.                     marked by C{location_of_fixed_flux} the value of C{u0} is
164           @type p0: scalar value on the domain (e.g. L{Data}).                     kept unchanged.
165           @param atol: absolute tolerance for the pressure          @type u0: vector value on the domain, e.g. L{Data}
166           @type atol: non-negative C{float}          @param p0: initial guess for the pressure. At locations in the domain
167           @param rtol: relative tolerance for the pressure                     marked by C{location_of_fixed_pressure} the value of C{p0}
168           @type rtol: non-negative C{float}                     is kept unchanged.
169           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}          @type p0: scalar value on the domain, e.g. L{Data}
170           @type sub_rtol: positive-negative C{float}          @param atol: absolute tolerance for the pressure
171           @param verbose: if set some information on iteration progress are printed          @type atol: non-negative C{float}
172           @type verbose: C{bool}          @param rtol: relative tolerance for the pressure
173           @param show_details:  if set information on the subiteration process are printed.          @type rtol: non-negative C{float}
174           @type show_details: C{bool}          @param sub_rtol: tolerance to be used in the sub iteration. It is
175           @return: flux and pressure                           recommended that M{sub_rtol<rtol*5.e-3}
176           @rtype: C{tuple} of L{Data}.          @type sub_rtol: positive-negative C{float}
177            @param verbose: if True information on iteration progress is printed
178           @note: The problem is solved as a least squares form          @type verbose: C{bool}
179            @param show_details: if True information on the sub-iteration process
180           M{(I+D^*D)u+Qp=D^*f+g}                               is printed
181           M{Q^*u+Q^*Qp=Q^*g}          @type show_details: C{bool}
182            @return: flux and pressure
183           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          @rtype: C{tuple} of L{Data}
184           We eliminate the flux form the problem by setting  
185            @note: the problem is solved in a least squares formulation:
186           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux  
187            M{(I+D^*D)u+Qp=D^*f+g}
188           form the first equation. Inserted into the second equation we get  
189            M{Q^*u+Q^*Qp=Q^*g}
190           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  
191                    where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the
192           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step          permeability M{k_{ij}}. We eliminate the flux from the problem by
193           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.          setting
194           """  
195           self.verbose=verbose          M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with M{u=u0} on C{location_of_fixed_flux}
196           self.show_details= show_details and self.verbose  
197           self.__pde_v.setTolerance(sub_rtol)          from the first equation. Inserted into the second equation we get
198           self.__pde_p.setTolerance(sub_rtol)  
199           u2=u0*self.__pde_v.getCoefficient("q")          M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with M{p=p0}
200           #          on C{location_of_fixed_pressure}
201           # first the reference velocity is calculated from  
202           #          which is solved using the PCG method (precondition is M{Q^*Q}).
203           #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)          In each iteration step PDEs with operator M{I+D^*D} and with M{Q^*Q}
204           #          need to be solved using a sub-iteration scheme.
205           self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)          """
206           u_ref=self.__pde_v.getSolution(verbose=show_details)          self.verbose=verbose
207           if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)          self.show_details= show_details and self.verbose
208           self.__pde_v.setValue(r=Data())          self.__pde_v.setTolerance(sub_rtol)
209           #          self.__pde_p.setTolerance(sub_rtol)
210           #   and then we calculate a reference pressure          u2=u0*self.__pde_v.getCoefficient("q")
211           #          #
212           #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)          # first the reference velocity is calculated from
213           #          #
214           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)          #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
215           p_ref=self.__pde_p.getSolution(verbose=self.show_details)          #
216           if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)          self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
217           self.__pde_p.setValue(r=Data())          u_ref=self.__pde_v.getSolution(verbose=show_details)
218           #          if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
219           #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref          self.__pde_v.setValue(r=Data())
220           #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref          #
221           #          #   and then we calculate a reference pressure
222           #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du          #
223           #          #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
224           #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref          #
225           #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref          self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
226           #          p_ref=self.__pde_p.getSolution(verbose=self.show_details)
227           #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)          if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
228           #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)          self.__pde_p.setValue(r=Data())
229           #          #
230           #  we use the (du,Qdp) to represent the resudual          #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
231           #  Q^*Q is a preconditioner          #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
232           #            #
233           #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>          #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
234           #          #
235           Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))          #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
236           norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))          #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
237           ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)          #
238           if not ATOL>0:          #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
239               raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, norm_rhs, atol)          #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
240           if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)          #
241           #          #  we use the (du,Qdp) to represent the resudual
242           #   caclulate the initial residual          #  Q^*Q is a preconditioner
243           #          #
244           self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())          #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
245           du=self.__pde_v.getSolution(verbose=show_details)          #
246           r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)          Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
247           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)          norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
248           util.saveVTK("d.vtu",p=dp,p_ref=p_ref)          ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
249           return u_ref+r[1],dp          if not ATOL>0:
250                        raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side = %e, atol =%e)."%(rtol, norm_rhs, atol)
251            if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
252            #
253            #   caclulate the initial residual
254            #
255            self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
256            du=self.__pde_v.getSolution(verbose=show_details)
257            r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
258            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)
259            util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
260            return u_ref+r[1],dp
261    
262      def __Aprod_PCG(self,p):      def __Aprod_PCG(self,p):
263            if self.show_details: print "DarcyFlux: Applying operator"          if self.show_details: print "DarcyFlux: Applying operator"
264            Qp=util.tensor_mult(self.__permeability,util.grad(p))          Qp=util.tensor_mult(self.__permeability,util.grad(p))
265            self.__pde_v.setValue(Y=Qp,X=Data())          self.__pde_v.setValue(Y=Qp,X=Data())
266            w=self.__pde_v.getSolution(verbose=self.show_details)          w=self.__pde_v.getSolution(verbose=self.show_details)
267            return ArithmeticTuple(-Qp,w)          return ArithmeticTuple(-Qp,w)
268    
269      def __inner_PCG(self,p,r):      def __inner_PCG(self,p,r):
270           a=util.tensor_mult(self.__permeability,util.grad(p))          a=util.tensor_mult(self.__permeability,util.grad(p))
271           out=-util.integrate(util.inner(a,r[0]+r[1]))          out=-util.integrate(util.inner(a,r[0]+r[1]))
272           return out          return out
273    
274      def __Msolve_PCG(self,r):      def __Msolve_PCG(self,r):
275            if self.show_details: print "DarcyFlux: Applying preconditioner"          if self.show_details: print "DarcyFlux: Applying preconditioner"
276            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]))
277            return self.__pde_p.getSolution(verbose=self.show_details)          return self.__pde_p.getSolution(verbose=self.show_details)
278    
279  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
280        """        """
281        solves        Represents and solves the problem
282    
283          M{-(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}}
284    
285            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}        M{u_{i,i}=0} and M{u=0} where C{fixed_u_mask}>0
                 u_{i,i}=0  
286    
287            u=0 where  fixed_u_mask>0        M{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  
288    
289        if surface_stress is not given 0 is assumed.        If C{surface_stress} is not given 0 is assumed.
290    
291        typical usage:        Typical usage::
292    
293              sp=StokesProblemCartesian(domain)            sp = StokesProblemCartesian(domain)
294              sp.setTolerance()            sp.setTolerance()
295              sp.initialize(...)            sp.initialize(...)
296              v,p=sp.solve(v0,p0)            v,p = sp.solve(v0,p0)
297        """        """
298        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
299           """           """
300           initialize the Stokes Problem           Initializes the Stokes Problem.
301    
302           @param domain: domain of the problem. The approximation order needs to be two.           @param domain: domain of the problem. The approximation order needs
303                            to be two.
304           @type domain: L{Domain}           @type domain: L{Domain}
305           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           @warning: The approximation order needs to be two otherwise you may
306                       see oscillations in the pressure.
307           """           """
308           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309           self.domain=domain           self.domain=domain
# Line 288  class StokesProblemCartesian(Homogeneous Line 312  class StokesProblemCartesian(Homogeneous
312           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
313           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
314           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
315                
316           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
317           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
318           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)           # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
# Line 300  class StokesProblemCartesian(Homogeneous Line 324  class StokesProblemCartesian(Homogeneous
324           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
325    
326        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()):
327          """           """
328          assigns values to the model parameters           Assigns values to the model parameters.
329    
330          @param f: external force           @param f: external force
331          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar           @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
332          @param fixed_u_mask: mask of locations with fixed velocity.           @param fixed_u_mask: mask of locations with fixed velocity
333          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar           @type fixed_u_mask: L{Vector} object on L{FunctionSpace}, L{Solution}
334          @param eta: viscosity                               or similar
335          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar           @param eta: viscosity
336          @param surface_stress: normal surface stress           @type eta: L{Scalar} object on L{FunctionSpace}, L{Function} or similar
337          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar           @param surface_stress: normal surface stress
338          @param stress: initial stress           @type surface_stress: L{Vector} object on L{FunctionSpace},
339      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar                                 L{FunctionOnBoundary} or similar
340          @note: All values needs to be set.           @param stress: initial stress
341             @type stress: L{Tensor} object on L{FunctionSpace}, L{Function} or
342          """                         similar
343          self.eta=eta           @note: All values need to be set.
344          A =self.__pde_u.createCoefficient("A")           """
345      self.__pde_u.setValue(A=Data())           self.eta=eta
346          for i in range(self.domain.getDim()):           A =self.__pde_u.createCoefficient("A")
347          for j in range(self.domain.getDim()):           self.__pde_u.setValue(A=Data())
348              A[i,j,j,i] += 1.           for i in range(self.domain.getDim()):
349              A[i,j,i,j] += 1.               for j in range(self.domain.getDim()):
350      self.__pde_prec.setValue(D=1/self.eta)                   A[i,j,j,i] += 1.
351          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)                   A[i,j,i,j] += 1.
352          self.__stress=stress           self.__pde_prec.setValue(D=1/self.eta)
353             self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
354             self.__stress=stress
355    
356        def B(self,v):        def B(self,v):
357          """           """
358          returns div(v)           Returns M{div(v)}.
359          @rtype: equal to the type of p           @return: M{div(v)}
360             @rtype: equal to the type of p
361    
362          @note: boundary conditions on p should be zero!           @note: Boundary conditions on p should be zero!
363          """           """
364          if self.show_details: print "apply divergence:"           if self.show_details: print "apply divergence:"
365          self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
366          self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setTolerance(self.getSubProblemTolerance())
367          return self.__pde_proj.getSolution(verbose=self.show_details)           return self.__pde_proj.getSolution(verbose=self.show_details)
368    
369        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
370           """           """
371           returns inner product of element p and Bv  (overwrite)           Returns inner product of element p and Bv (overwrite).
372            
373           @type p: equal to the type of p           @type p: equal to the type of p
374           @type Bv: equal to the type of result of operator B           @type Bv: equal to the type of result of operator B
375           @rtype: C{float}           @return: inner product of p and Bv
   
376           @rtype: equal to the type of p           @rtype: equal to the type of p
377           """           """
378           s0=util.interpolate(p,Function(self.domain))           s0=util.interpolate(p,Function(self.domain))
# Line 355  class StokesProblemCartesian(Homogeneous Line 381  class StokesProblemCartesian(Homogeneous
381    
382        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
383           """           """
384           returns inner product of element p0 and p1  (overwrite)           Returns inner product of element p0 and p1 (overwrite).
385            
386           @type p0: equal to the type of p           @type p0: equal to the type of p
387           @type p1: equal to the type of p           @type p1: equal to the type of p
388           @rtype: C{float}           @return: inner product of p0 and p1
   
389           @rtype: equal to the type of p           @rtype: equal to the type of p
390           """           """
391           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0/self.eta,Function(self.domain))
# Line 369  class StokesProblemCartesian(Homogeneous Line 394  class StokesProblemCartesian(Homogeneous
394    
395        def inner_v(self,v0,v1):        def inner_v(self,v0,v1):
396           """           """
397           returns inner product of two element v0 and v1  (overwrite)           Returns inner product of two elements v0 and v1 (overwrite).
398            
399           @type v0: equal to the type of v           @type v0: equal to the type of v
400           @type v1: equal to the type of v           @type v1: equal to the type of v
401           @rtype: C{float}           @return: inner product of v0 and v1
   
402           @rtype: equal to the type of v           @rtype: equal to the type of v
403           """           """
404       gv0=util.grad(v0)           gv0=util.grad(v0)
405       gv1=util.grad(v1)           gv1=util.grad(v1)
406           return util.integrate(util.inner(gv0,gv1))           return util.integrate(util.inner(gv0,gv1))
407    
408        def solve_A(self,u,p):        def solve_A(self,u,p):
409           """           """
410           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           Solves M{Av=f-Au-B^*p} (v=0 on fixed_u_mask).
411           """           """
412           if self.show_details: print "solve for velocity:"           if self.show_details: print "solve for velocity:"
413           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
# Line 392  class StokesProblemCartesian(Homogeneous Line 416  class StokesProblemCartesian(Homogeneous
416           else:           else:
417              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))
418           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_u.getSolution(verbose=self.show_details)
419           return  out           return out
420    
421        def solve_prec(self,p):        def solve_prec(self,p):
422             """
423             Applies the preconditioner.
424             """
425           if self.show_details: print "apply preconditioner:"           if self.show_details: print "apply preconditioner:"
426           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
427           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
428           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
429           return q           return q
430    

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

  ViewVC Help
Powered by ViewVC 1.1.26