/[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 2445 by gross, Fri May 29 03:23:25 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
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        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):
294             return util.integrate(util.inner(self.__Q(p), r))
295    
296        def __Msolve_PCG(self,r):
297              self.__pde_p.setTolerance(self.getSubProblemTolerance())
298              if self.show_details: print "DarcyFlux: Applying preconditioner"
299              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300              # 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            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
335                  u_{i,i}=0                  u_{i,i}=0
336    
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=surface_stress            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 give 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
352    
353             @param domain: domain of the problem. The approximation order needs to be two.
354             @type domain: L{Domain}
355             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
356             """
357           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
358           self.domain=domain           self.domain=domain
359           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
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)
363           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # 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)
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()):  
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
379    
380            @param f: external force
381            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
382            @param fixed_u_mask: mask of locations with fixed velocity.
383            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
384            @param eta: viscosity
385            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
386            @param surface_stress: normal surface stress
387            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
388            @param stress: initial stress
389        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
390            @note: All values needs to be set.
391    
392            """
393          self.eta=eta          self.eta=eta
394          A =self.__pde_u.createCoefficientOfGeneralPDE("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
405    
406         def Bv(self,v):
407             """
408             returns inner product of element p and div(v)
409    
410             @param p: a pressure increment
411             @param v: a residual
412             @return: inner product of element p and div(v)
413             @rtype: C{float}
414             """
415             self.__pde_proj.setValue(Y=-util.div(v))
416             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
417             return self.__pde_proj.getSolution()
418    
419         def inner_pBv(self,p,Bv):
420             """
421             returns inner product of element p and Bv=-div(v)
422    
423             @param p: a pressure increment
424             @param v: a residual
425             @return: inner product of element p and Bv=-div(v)
426             @rtype: C{float}
427             """
428             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
429    
430         def inner_p(self,p0,p1):
431             """
432             Returns inner product of p0 and p1
433    
434        def B(self,arg):           @param p0: a pressure
435           d=util.div(arg)           @param p1: a pressure
436           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
437           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
438           return self.__pde_proj.getSolution(verbose=self.show_details)           """
439             s0=util.interpolate(p0/self.eta,Function(self.domain))
440        def inner(self,p0,p1):           s1=util.interpolate(p1/self.eta,Function(self.domain))
          s0=util.interpolate(p0,Function(self.domain))  
          s1=util.interpolate(p1,Function(self.domain))  
441           return util.integrate(s0*s1)           return util.integrate(s0*s1)
442    
443        def inner_a(self,a0,a1):       def norm_v(self,v):
444           p0=util.interpolate(a0[1],Function(self.domain))           """
445           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
446           alfa=(1/self.vol)*util.integrate(p0)  
447           beta=(1/self.vol)*util.integrate(p1)           @param v: a velovity
448       v0=util.grad(a0[0])           @return: norm of v
449       v1=util.grad(a1[0])           @rtype: non-negative C{float}
450           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           """
451             return util.sqrt(util.integrate(util.length(util.grad(v))))
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
452    
453        def solve_A(self,u,p):       def getV(self, p, v0):
454           """           """
455           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p (overwrite)
456    
457             @param p: a pressure
458             @param v0: a initial guess for the value v to return.
459             @return: v given as M{v= A^{-1} (f-B^*p)}
460           """           """
461           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
462           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)
463           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
464                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
465             else:
466                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
467             out=self.__pde_u.getSolution(verbose=self.show_details)
468             return  out
469    
470         def norm_Bv(self,Bv):
471            """
472            Returns Bv (overwrite).
473    
474            @rtype: equal to the type of p
475            @note: boundary conditions on p should be zero!
476            """
477            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
478    
479         def solve_AinvBt(self,p):
480             """
481             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
482    
483        def solve_prec(self,p):           @param p: a pressure increment
484       #proj=Projector(domain=self.domain, reduce = True, fast=False)           @return: the solution of M{Av=B^*p}
485           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           @note: boundary conditions on v should be zero!
486           self.__pde_prec.setValue(Y=p)           """
487           q=self.__pde_prec.getSolution(verbose=self.show_details)           self.__pde_u.setTolerance(self.getSubProblemTolerance())
488           return q           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
489             out=self.__pde_u.getSolution(verbose=self.show_details)
490             return  out
491    
492        def solve_prec1(self,p):       def solve_prec(self,Bv):
493       #proj=Projector(domain=self.domain, reduce = True, fast=False)           """
494           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
495           self.__pde_prec.setValue(Y=p)           with accuracy L{self.getSubProblemTolerance()}
          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  
496    
497             @param v: velocity increment
498             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
499             @note: boundary conditions on p are zero.
500             """
501             self.__pde_prec.setValue(Y=Bv)
502             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
503             return self.__pde_prec.getSolution(verbose=self.show_details)
504    
505    # rename solve_prec and change argument v to Bv
506    # chnage the argument of inner_pBv to v->Bv
507    # add Bv
508    # inner p still needed?
509    # change norm_Bv argument to Bv

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

  ViewVC Help
Powered by ViewVC 1.1.26