/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1673 by gross, Thu Jul 24 22:28:50 2008 UTC revision 2354 by gross, Wed Apr 1 04:01:58 2009 UTC
# Line 1  Line 1 
1  # $Id:$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2008 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8  #  # Licensed under the Open Software License version 3.0
9  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
10  #  #
11    ########################################################
12    
13    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
16    Primary Business: Queensland, Australia"""
17    __license__="""Licensed under the Open Software License version 3.0
18    http://www.opensource.org/licenses/osl-3.0.php"""
19    __url__="https://launchpad.net/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39  class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  class DarcyFlow(object):
40        """      """
41        solves      solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        M{u_{i,i} = f}
45    
46        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.
49        """
50    
51        def __init__(self, domain, weight=None, useReduced=False):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56            """
57            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)
64            if useReduced: self.__pde_v.setReducedOrderOn()
65            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)
69            self.__pde_p.setSymmetryOn()
70            if useReduced: self.__pde_p.setReducedOrderOn()
71            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72            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):
78            """
79            assigns values to model parameters
80    
81            @param f: volumetic sources/sinks
82            @type f: scalar value on the domain (e.g. L{Data})
83            @param g: flux sources/sinks
84            @type g: vector values on the domain (e.g. L{Data})
85            @param location_of_fixed_pressure: mask for locations where pressure is fixed
86            @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.
88            @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
90                                 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.
92            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93    
94            @note: the values of parameters which are not set by calling C{setValue} are not altered.
95            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96                   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.
98            """
99            if f !=None:
100               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101               if f.isEmpty():
102                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103               else:
104                   if f.getRank()>0: raise ValueError,"illegal rank of f."
105               self.__f=f
106            if g !=None:
107               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108               if g.isEmpty():
109                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110               else:
111                 if not g.getShape()==(self.domain.getDim(),):
112                   raise ValueError,"illegal shape of g"
113               self.__g=g
114    
115            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117    
118            if permeability!=None:
119               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120               if perm.getRank()==0:
121                   perm=perm*util.kronecker(self.domain.getDim())
122               elif perm.getRank()==1:
123                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125               elif perm.getRank()==2:
126                  pass
127               else:
128                  raise ValueError,"illegal rank of permeability."
129               self.__permeability=perm
130               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            @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142            """
143            if rtol<0:
144                raise ValueError,"Relative tolerance needs to be non-negative."
145            self.__rtol=rtol
146        def getTolerance(self):
147            """
148            returns the relative tolerance
149    
150            @return: current relative tolerance
151            @rtype: C{float}
152            """
153            return self.__rtol
154    
155        def setAbsoluteTolerance(self,atol=0.):
156            """
157            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158    
159            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            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i      def setSubProblemTolerance(self,rtol=None):
179                  u_{i,i}=0           """
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            u=0 where  fixed_u_mask>0           @param rtol: relative tolerence
184            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress           @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        if surface_stress is not give 0 is assumed.      def getSubProblemTolerance(self):
196             """
197             Returns the subproblem reduction factor.
198    
199        typical usage:           @return: subproblem reduction factor
200             @rtype: C{float}
201             """
202             return self.__sub_tol
203    
204              sp=StokesProblemCartesian(domain)      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205              sp.setTolerance()           """
206              sp.initialize(...)           solves the problem.
             v,p=sp.solve(v0,p0)  
       """  
       def __init__(self,domain,**kwargs):  
          HomogeneousSaddlePointProblem.__init__(self,**kwargs)  
          self.domain=domain  
          self.vol=util.integrate(1.,Function(self.domain))  
          self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())  
          self.__pde_u.setSymmetryOn()  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
          # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)  
          # self.__pde_proj.setReducedOrderOn()  
          # self.__pde_proj.setSymmetryOn()  
          # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)  
207    
208        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):           The iteration is terminated if the residual norm is less then self.getTolerance().
         self.eta=eta  
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
         # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
209    
210          # self.__pde_proj.setValue(D=1/eta)           @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          # self.__pde_proj.setValue(Y=1.)           @type u0: vector value on the domain (e.g. L{Data}).
212          # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))           @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          self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))           @type p0: scalar value on the domain (e.g. L{Data}).
214             @param verbose: if set some information on iteration progress are printed
215        def B(self,arg):           @type verbose: C{bool}
216           a=util.div(arg, ReducedFunction(self.domain))           @param show_details:  if set information on the subiteration process are printed.
217           return a-util.integrate(a)/self.vol           @type show_details: C{bool}
218             @return: flux and pressure
219        def inner(self,p0,p1):           @rtype: C{tuple} of L{Data}.
          return util.integrate(p0*p1)  
          # return util.integrate(1/self.__inv_eta**2*p0*p1)  
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
220    
221        def solve_A(self,u,p):           @note: The problem is solved as a least squares form
222           """  
223           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           M{(I+D^*D)u+Qp=D^*f+g}
224           """           M{Q^*u+Q^*Qp=Q^*g}
225           self.__pde_u.setTolerance(self.getSubProblemTolerance())  
226           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227           return  self.__pde_u.getSolution(verbose=self.show_details)           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
230    
231        def solve_prec(self,p):           form the first equation. Inserted into the second equation we get
         a=self.__inv_eta*p  
         return a-util.integrate(a)/self.vol  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
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
234    
235             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.
237             """
238             self.verbose=verbose or True
239             self.show_details= show_details and self.verbose
240             rtol=self.getTolerance()
241             atol=self.getAbsoluteTolerance()
242             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243    
244             num_corrections=0
245             converged=False
246             p=p0
247             norm_r=None
248             while not converged:
249                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250                   Qp=self.__Q(p)
251                   norm_v=self.__L2(v)
252                   norm_Qp=self.__L2(Qp)
253                   if norm_v == 0.:
254                      if norm_Qp == 0.:
255                         return v,p
256                      else:
257                        fac=norm_Qp
258                   else:
259                      if norm_Qp == 0.:
260                        fac=norm_v
261                      else:
262                        fac=2./(1./norm_v+1./norm_Qp)
263                   ATOL=(atol+rtol*fac)
264                   if self.verbose:
265                        print "DarcyFlux: L2 norm of v = %e."%norm_v
266                        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    #
277    #              
278    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
279    #               #===========================================================================
280    #               norm_r_hat=self.__L2(r_hat)
281    #               norm_v=self.__L2(v)
282    #               norm_g=self.__L2(g)
283    #               norm_gv=self.__L2(g-v)
284    #               norm_Qp=self.__L2(Qp)
285    #               norm_gQp=self.__L2(g-Qp)
286    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
287    #               fac=min(norm_v,norm_Qp,norm_gv)
288    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
289    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
290    #               if r!=None:
291    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
292    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
293    #                   self.setSubProblemTolerance(sub_tol)
294    #                   print "subtol_new=",self.getSubProblemTolerance()
295    #               print "norm_v = ",norm_v
296    #               print "norm_gv = ",norm_gv
297    #               print "norm_Qp = ",norm_Qp
298    #               print "norm_gQp = ",norm_gQp
299    #               print "norm_g = ",norm_g
300    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
301    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
302    #               if fac == 0:
303    #                   if self.verbose: print "DarcyFlux: trivial case!"
304    #                   return v,p
305    #               #===============================================================================
306    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
307    #               # norm_Qp=self.__L2(Qp)
308    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
309    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
310    #
311    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
312    #               if self.verbose:
313    #                   print "DarcyFlux: residual = %e"%norm_r_hat
314    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
315    #               if norm_r_hat <= ATOL:
316    #                   print "DarcyFlux: iteration finalized."
317    #                   converged=True
318    #               else:
319    #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)
320    #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)
321    #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
322    #               print "norm_r =",norm_r
323    #         return v,p
324        def __L2(self,v):
325             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
326    
327        def __Q(self,p):
328              return util.tensor_mult(self.__permeability,util.grad(p))
329    
330        def __Aprod(self,dp):
331              self.__pde_v.setTolerance(self.getSubProblemTolerance())
332              if self.show_details: print "DarcyFlux: Applying operator"
333              Qdp=self.__Q(dp)
334              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
335              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
336              return Qdp+du
337        def __inner_GMRES(self,r,s):
338             return util.integrate(util.inner(r,s))
339    
340        def __inner_PCG(self,p,r):
341             return util.integrate(util.inner(self.__Q(p), r))
342    
343        def __Msolve_PCG(self,r):
344              self.__pde_p.setTolerance(self.getSubProblemTolerance())
345              if self.show_details: print "DarcyFlux: Applying preconditioner"
346              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
347              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
348    
349        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
350            """
351            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
352            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
353            Note that C{g} and C{f} are used, see L{setValue}.
354    
355            @param p: pressure.
356            @type p: scalar value on the domain (e.g. L{Data}).
357            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
358            @type fixed_flux: vector values on the domain (e.g. L{Data}).
359            @param tol: relative tolerance to be used.
360            @type tol: positive C{float}.
361            @return: flux
362            @rtype: L{Data}
363            @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}}
364                   for the permeability M{k_{ij}}
365            """
366            self.__pde_v.setTolerance(self.getSubProblemTolerance())
367            g=self.__g
368            f=self.__f
369            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
370            if p == None:
371               self.__pde_v.setValue(Y=g)
372            else:
373               self.__pde_v.setValue(Y=g-self.__Q(p))
374            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
375    
376  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
377        """       """
378        solves       solves
379    
380            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
381                  u_{i,i}=0                  u_{i,i}=0
382    
383            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
384            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
385    
386        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
387    
388        typical usage:       typical usage:
389    
390              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
391              sp.setTolerance()              sp.setTolerance()
392              sp.initialize(...)              sp.initialize(...)
393              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
394        """       """
395        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
396             """
397             initialize the Stokes Problem
398    
399             @param domain: domain of the problem. The approximation order needs to be two.
400             @type domain: L{Domain}
401             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
402             """
403           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
404           self.domain=domain           self.domain=domain
405           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
406           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())
407           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
408             # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
409           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
410                
411           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
412           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
413             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
414           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
415    
416           self.__pde_proj=LinearPDE(domain)       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
417           self.__pde_proj.setReducedOrderOn()          """
418           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
419           self.__pde_proj.setValue(D=1.)  
420            @param f: external force
421            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
422            @param fixed_u_mask: mask of locations with fixed velocity.
423            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
424            @param eta: viscosity
425            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
426            @param surface_stress: normal surface stress
427            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
428            @param stress: initial stress
429        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
430            @note: All values needs to be set.
431    
432        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          """
433          self.eta=eta          self.eta=eta
434          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
435      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
436          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
437          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
438              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
439              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
440      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
441          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)
442            self.__f=f
443            self.__surface_stress=surface_stress
444            self.__stress=stress
445    
446        def B(self,arg):       def inner_pBv(self,p,v):
447           d=util.div(arg)           """
448           self.__pde_proj.setValue(Y=d)           returns inner product of element p and div(v)
449           self.__pde_proj.setTolerance(self.getSubProblemTolerance())  
450           return self.__pde_proj.getSolution(verbose=self.show_details)           @param p: a pressure increment
451             @param v: a residual
452        def inner(self,p0,p1):           @return: inner product of element p and div(v)
453           s0=util.interpolate(p0,Function(self.domain))           @rtype: C{float}
454           s1=util.interpolate(p1,Function(self.domain))           """
455             return util.integrate(-p*util.div(v))
456    
457         def inner_p(self,p0,p1):
458             """
459             Returns inner product of p0 and p1
460    
461             @param p0: a pressure
462             @param p1: a pressure
463             @return: inner product of p0 and p1
464             @rtype: C{float}
465             """
466             s0=util.interpolate(p0/self.eta,Function(self.domain))
467             s1=util.interpolate(p1/self.eta,Function(self.domain))
468           return util.integrate(s0*s1)           return util.integrate(s0*s1)
469    
470        def inner_a(self,a0,a1):       def norm_v(self,v):
471           p0=util.interpolate(a0[1],Function(self.domain))           """
472           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
          alfa=(1/self.vol)*util.integrate(p0)  
          beta=(1/self.vol)*util.integrate(p1)  
      v0=util.grad(a0[0])  
      v1=util.grad(a1[0])  
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
   
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
473    
474        def solve_A(self,u,p):           @param v: a velovity
475             @return: norm of v
476             @rtype: non-negative C{float}
477           """           """
478           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return util.sqrt(util.integrate(util.length(util.grad(v))))
479    
480         def getV(self, p, v0):
481             """
482             return the value for v for a given p (overwrite)
483    
484             @param p: a pressure
485             @param v0: a initial guess for the value v to return.
486             @return: v given as M{v= A^{-1} (f-B^*p)}
487           """           """
488           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
489           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
490           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
491                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
492             else:
493                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
494             out=self.__pde_u.getSolution(verbose=self.show_details)
495             return  out
496    
497    
498        def solve_prec(self,p):           raise NotImplementedError,"no v calculation implemented."
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
499    
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
      q0=util.interpolate(q,Function(self.domain))  
          print util.inf(q*q0),util.sup(q*q0)  
          q-=(1/self.vol)*util.integrate(q0)  
          print util.inf(q*q0),util.sup(q*q0)  
          return q  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
500    
501         def norm_Bv(self,v):
502            """
503            Returns Bv (overwrite).
504    
505            @rtype: equal to the type of p
506            @note: boundary conditions on p should be zero!
507            """
508            return util.sqrt(util.integrate(util.div(v)**2))
509    
510         def solve_AinvBt(self,p):
511             """
512             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
513    
514             @param p: a pressure increment
515             @return: the solution of M{Av=B^*p}
516             @note: boundary conditions on v should be zero!
517             """
518             self.__pde_u.setTolerance(self.getSubProblemTolerance())
519             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
520             out=self.__pde_u.getSolution(verbose=self.show_details)
521             return  out
522    
523         def solve_precB(self,v):
524             """
525             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
526             with accuracy L{self.getSubProblemTolerance()} (overwrite).
527    
528             @param v: velocity increment
529             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
530             @note: boundary conditions on p are zero.
531             """
532             self.__pde_prec.setValue(Y=-util.div(v))
533             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
534             return self.__pde_prec.getSolution(verbose=self.show_details)

Legend:
Removed from v.1673  
changed lines
  Added in v.2354

  ViewVC Help
Powered by ViewVC 1.1.26