/[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 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1 
1  # $Id:$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2010 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-2010 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
23    
24  @var __author__: name of author  :var __author__: name of author
25  @var __copyright__: copyrights  :var __copyright__: copyrights
26  @var __license__: licence agreement  :var __license__: licence agreement
27  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
28  @var __version__: version  :var __version__: version
29  @var __date__: date of the version  :var __date__: date of the version
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, SolverOptions
37  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38    
39    class DarcyFlow(object):
40        """
41        solves the problem
42    
43        *u_i+k_{ij}*p_{,j} = g_i*
44        *u_{i,i} = f*
45    
46        where *p* represents the pressure and *u* the Darcy flux. *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, adaptSubTolerance=True):
52            """
53            initializes the Darcy flux problem
54            :param domain: domain of the problem
55            :type domain: `Domain`
56        :param useReduced: uses reduced oreder on flux and pressure
57        :type useReduced: ``bool``
58        :param adaptSubTolerance: switches on automatic subtolerance selection
59        :type adaptSubTolerance: ``bool``  
60            """
61            self.domain=domain
62            if weight == None:
63               s=self.domain.getSize()
64               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65               # self.__l=(3.*util.longestEdge(self.domain))**2
66               #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
67            else:
68               self.__l=weight
69            self.__pde_v=LinearPDESystem(domain)
70            if useReduced: self.__pde_v.setReducedOrderOn()
71            self.__pde_v.setSymmetryOn()
72            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
73            self.__pde_p=LinearSinglePDE(domain)
74            self.__pde_p.setSymmetryOn()
75            if useReduced: self.__pde_p.setReducedOrderOn()
76            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
77            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78            self.setTolerance()
79            self.setAbsoluteTolerance()
80        self.__adaptSubTolerance=adaptSubTolerance
81        self.verbose=False
82        def getSolverOptionsFlux(self):
83        """
84        Returns the solver options used to solve the flux problems
85        
86        *(I+D^*D)u=F*
87        
88        :return: `SolverOptions`
89        """
90        return self.__pde_v.getSolverOptions()
91        def setSolverOptionsFlux(self, options=None):
92        """
93        Sets the solver options used to solve the flux problems
94        
95        *(I+D^*D)u=F*
96        
97        If ``options`` is not present, the options are reset to default
98        :param options: `SolverOptions`
99        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
100        """
101        return self.__pde_v.setSolverOptions(options)
102        def getSolverOptionsPressure(self):
103        """
104        Returns the solver options used to solve the pressure problems
105        
106        *(Q^*Q)p=Q^*G*
107        
108        :return: `SolverOptions`
109        """
110        return self.__pde_p.getSolverOptions()
111        def setSolverOptionsPressure(self, options=None):
112        """
113        Sets the solver options used to solve the pressure problems
114        
115        *(Q^*Q)p=Q^*G*
116        
117        If ``options`` is not present, the options are reset to default
118        :param options: `SolverOptions`
119        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
120        """
121        return self.__pde_p.setSolverOptions(options)
122    
123        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124            """
125            assigns values to model parameters
126    
127            :param f: volumetic sources/sinks
128            :type f: scalar value on the domain (e.g. `Data`)
129            :param g: flux sources/sinks
130            :type g: vector values on the domain (e.g. `Data`)
131            :param location_of_fixed_pressure: mask for locations where pressure is fixed
132            :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
133            :param location_of_fixed_flux:  mask for locations where flux is fixed.
134            :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
135            :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
136                                 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
137                                 ``v`` on the main diagonal is used.
138            :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
139    
140            :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141            :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
142                   or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
143                   is along the *x_i* axis.
144            """
145            if f !=None:
146               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
147               if f.isEmpty():
148                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
149               else:
150                   if f.getRank()>0: raise ValueError,"illegal rank of f."
151               self.__f=f
152            if g !=None:
153               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
154               if g.isEmpty():
155                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
156               else:
157                 if not g.getShape()==(self.domain.getDim(),):
158                   raise ValueError,"illegal shape of g"
159               self.__g=g
160    
161            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
162            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
163    
164            if permeability!=None:
165               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
166               if perm.getRank()==0:
167                   perm=perm*util.kronecker(self.domain.getDim())
168               elif perm.getRank()==1:
169                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
170                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
171               elif perm.getRank()==2:
172                  pass
173               else:
174                  raise ValueError,"illegal rank of permeability."
175               self.__permeability=perm
176               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
177    
178        def setTolerance(self,rtol=1e-4):
179            """
180            sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
181    
182            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183    
184            where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
185    
186            :param rtol: relative tolerance for the pressure
187            :type rtol: non-negative ``float``
188            """
189            if rtol<0:
190                raise ValueError,"Relative tolerance needs to be non-negative."
191            self.__rtol=rtol
192        def getTolerance(self):
193            """
194            returns the relative tolerance
195    
196            :return: current relative tolerance
197            :rtype: ``float``
198            """
199            return self.__rtol
200    
201        def setAbsoluteTolerance(self,atol=0.):
202            """
203            sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204    
205            *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206    
207            where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
208    
209            :param atol: absolute tolerance for the pressure
210            :type atol: non-negative ``float``
211            """
212            if atol<0:
213                raise ValueError,"Absolute tolerance needs to be non-negative."
214            self.__atol=atol
215        def getAbsoluteTolerance(self):
216           """
217           returns the absolute tolerance
218          
219           :return: current absolute tolerance
220           :rtype: ``float``
221           """
222           return self.__atol
223        def getSubProblemTolerance(self):
224        """
225        Returns a suitable subtolerance
226        @type: ``float``
227        """
228        return max(util.EPSILON**(0.75),self.getTolerance()**2)
229        def setSubProblemTolerance(self):
230             """
231             Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
232             """
233         if self.__adaptSubTolerance:
234             sub_tol=self.getSubProblemTolerance()
235                 self.getSolverOptionsFlux().setTolerance(sub_tol)
236             self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
237             self.getSolverOptionsPressure().setTolerance(sub_tol)
238             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
239             if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
240    
241        def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
242             """
243             solves the problem.
244    
245             The iteration is terminated if the residual norm is less then self.getTolerance().
246    
247             :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
248             :type u0: vector value on the domain (e.g. `Data`).
249             :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
250             :type p0: scalar value on the domain (e.g. `Data`).
251             :param verbose: if set some information on iteration progress are printed
252             :type verbose: ``bool``
253             :return: flux and pressure
254             :rtype: ``tuple`` of `Data`.
255    
256             :note: The problem is solved as a least squares form
257    
258             *(I+D^*D)u+Qp=D^*f+g*
259             *Q^*u+Q^*Qp=Q^*g*
260    
261             where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
262             We eliminate the flux form the problem by setting
263    
264             *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux
265    
266             form the first equation. Inserted into the second equation we get
267    
268             *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure
269    
270             which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
271             PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
272             """
273             self.verbose=verbose
274             rtol=self.getTolerance()
275             atol=self.getAbsoluteTolerance()
276         self.setSubProblemTolerance()
277             num_corrections=0
278             converged=False
279             p=p0
280             norm_r=None
281             while not converged:
282                   v=self.getFlux(p, fixed_flux=u0)
283                   Qp=self.__Q(p)
284                   norm_v=self.__L2(v)
285                   norm_Qp=self.__L2(Qp)
286                   if norm_v == 0.:
287                      if norm_Qp == 0.:
288                         return v,p
289                      else:
290                        fac=norm_Qp
291                   else:
292                      if norm_Qp == 0.:
293                        fac=norm_v
294                      else:
295                        fac=2./(1./norm_v+1./norm_Qp)
296                   ATOL=(atol+rtol*fac)
297                   if self.verbose:
298                        print "DarcyFlux: L2 norm of v = %e."%norm_v
299                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
300                        print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
301                        print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
302                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
303                   if norm_r == None or norm_r>ATOL:
304                       if num_corrections>max_num_corrections:
305                             raise ValueError,"maximum number of correction steps reached."
306                       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)
307                       num_corrections+=1
308                   else:
309                       converged=True
310             return v,p
311        def __L2(self,v):
312             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
313    
314        def __Q(self,p):
315              return util.tensor_mult(self.__permeability,util.grad(p))
316    
317        def __Aprod(self,dp):
318              if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
319              Qdp=self.__Q(dp)
320              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
321              du=self.__pde_v.getSolution()
322              # self.__pde_v.getOperator().saveMM("proj.mm")
323              return Qdp+du
324        def __inner_GMRES(self,r,s):
325             return util.integrate(util.inner(r,s))
326    
327        def __inner_PCG(self,p,r):
328             return util.integrate(util.inner(self.__Q(p), r))
329    
330        def __Msolve_PCG(self,r):
331          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
332              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
333              # self.__pde_p.getOperator().saveMM("prec.mm")
334              return self.__pde_p.getSolution()
335    
336        def getFlux(self,p=None, fixed_flux=Data()):
337            """
338            returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
339            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
340            Note that ``g`` and ``f`` are used, see `setValue`.
341    
342            :param p: pressure.
343            :type p: scalar value on the domain (e.g. `Data`).
344            :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
345            :type fixed_flux: vector values on the domain (e.g. `Data`).
346            :return: flux
347            :rtype: `Data`
348            :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
349                   for the permeability *k_{ij}*
350            """
351        self.setSubProblemTolerance()
352            g=self.__g
353            f=self.__f
354            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
355            if p == None:
356               self.__pde_v.setValue(Y=g)
357            else:
358               self.__pde_v.setValue(Y=g-self.__Q(p))
359            return self.__pde_v.getSolution()
360    
361  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
362        """       """
363        solves       solves
364    
365            -(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}
366                  u_{i,i}=0                  u_{i,i}=0
367    
368            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
369            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
370    
371        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
372    
373        typical usage:       typical usage:
374    
375              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
376              sp.setTolerance()              sp.setTolerance()
377              sp.initialize(...)              sp.initialize(...)
378              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
379        """       """
380        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
381             """
382             initialize the Stokes Problem
383    
384             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
385             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
386             with macro elements for the pressure.
387    
388             :param domain: domain of the problem.
389             :type domain: `Domain`
390             """
391           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
392           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
393           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())
394           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
395           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)      
               
396           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
397           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
398           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
399    
400           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
401           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
402         self.__pde_proj.setValue(D=1)
403           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
404    
405        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
406          self.eta=eta           """
407          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
408      self.__pde_u.setValue(A=Data())      
409          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
410          for j in range(self.domain.getDim()):       """
411              A[i,j,j,i] += 1.       return self.__pde_u.getSolverOptions()
412              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
413      self.__pde_prec.setValue(D=1/self.eta)           """
414          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)       set the solver options for solving the equation for velocity.
415        
416        def B(self,arg):       :param options: new solver  options
417           d=util.div(arg)       :type options: `SolverOptions`
418           self.__pde_proj.setValue(Y=d)       """
419           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setSolverOptions(options)
420           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
421             """
422         returns the solver options used  solve the equation for pressure.
423         :rtype: `SolverOptions`
424         """
425         return self.__pde_prec.getSolverOptions()
426         def setSolverOptionsPressure(self, options=None):
427             """
428         set the solver options for solving the equation for pressure.
429         :param options: new solver  options
430         :type options: `SolverOptions`
431         """
432         self.__pde_prec.setSolverOptions(options)
433    
434        def inner(self,p0,p1):       def setSolverOptionsDiv(self, options=None):
435             """
436         set the solver options for solving the equation to project the divergence of
437         the velocity onto the function space of presure.
438        
439         :param options: new solver options
440         :type options: `SolverOptions`
441         """
442         self.__pde_proj.setSolverOptions(options)
443         def getSolverOptionsDiv(self):
444             """
445         returns the solver options for solving the equation to project the divergence of
446         the velocity onto the function space of presure.
447        
448         :rtype: `SolverOptions`
449         """
450         return self.__pde_proj.getSolverOptions()
451    
452         def updateStokesEquation(self, v, p):
453             """
454             updates the Stokes equation to consider dependencies from ``v`` and ``p``
455             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
456             """
457             pass
458         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
459            """
460            assigns new values to the model parameters.
461    
462            :param f: external force
463            :type f: `Vector` object in `FunctionSpace` `Function` or similar
464            :param fixed_u_mask: mask of locations with fixed velocity.
465            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
466            :param eta: viscosity
467            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
468            :param surface_stress: normal surface stress
469            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
470            :param stress: initial stress
471        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
472            """
473            if eta !=None:
474                k=util.kronecker(self.domain.getDim())
475                kk=util.outer(k,k)
476                self.eta=util.interpolate(eta, Function(self.domain))
477            self.__pde_prec.setValue(D=1/self.eta)
478                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
479            if restoration_factor!=None:
480                n=self.domain.getNormal()
481                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
482            if fixed_u_mask!=None:
483                self.__pde_u.setValue(q=fixed_u_mask)
484            if f!=None: self.__f=f
485            if surface_stress!=None: self.__surface_stress=surface_stress
486            if stress!=None: self.__stress=stress
487    
488         def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
489            """
490            assigns values to the model parameters
491    
492            :param f: external force
493            :type f: `Vector` object in `FunctionSpace` `Function` or similar
494            :param fixed_u_mask: mask of locations with fixed velocity.
495            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
496            :param eta: viscosity
497            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
498            :param surface_stress: normal surface stress
499            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
500            :param stress: initial stress
501        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
502            """
503            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
504    
505         def Bv(self,v,tol):
506             """
507             returns inner product of element p and div(v)
508    
509             :param v: a residual
510             :return: inner product of element p and div(v)
511             :rtype: ``float``
512             """
513             self.__pde_proj.setValue(Y=-util.div(v))
514         self.getSolverOptionsDiv().setTolerance(tol)
515         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
516             out=self.__pde_proj.getSolution()
517             return out
518    
519         def inner_pBv(self,p,Bv):
520             """
521             returns inner product of element p and Bv=-div(v)
522    
523             :param p: a pressure increment
524             :param Bv: a residual
525             :return: inner product of element p and Bv=-div(v)
526             :rtype: ``float``
527             """
528             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
529    
530         def inner_p(self,p0,p1):
531             """
532             Returns inner product of p0 and p1
533    
534             :param p0: a pressure
535             :param p1: a pressure
536             :return: inner product of p0 and p1
537             :rtype: ``float``
538             """
539           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
540           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
541           return util.integrate(s0*s1)           return util.integrate(s0*s1)
542    
543        def inner_a(self,a0,a1):       def norm_v(self,v):
544           p0=util.interpolate(a0[1],Function(self.domain))           """
545           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
546           alfa=(1/self.vol)*util.integrate(p0)  
547           beta=(1/self.vol)*util.integrate(p1)           :param v: a velovity
548       v0=util.grad(a0[0])           :return: norm of v
549       v1=util.grad(a1[0])           :rtype: non-negative ``float``
550       #print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))           """
551           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
552    
553    
554        def getStress(self,u):       def getDV(self, p, v, tol):
555           mg=util.grad(u)           """
556           return 2.*self.eta*util.symmetric(mg)           return the value for v for a given p (overwrite)
557    
558        def solve_A(self,u,p):           :param p: a pressure
559           """           :param v: a initial guess for the value v to return.
560           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :return: dv given as *Adv=(f-Av-B^*p)*
561           """           """
562           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.updateStokesEquation(v,p)
563           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
564           return  self.__pde_u.getSolution(verbose=self.show_details)       self.getSolverOptionsVelocity().setTolerance(tol)
565         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
566             if self.__stress.isEmpty():
567        def solve_prec(self,p):              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
568       #proj=Projector(domain=self.domain, reduce = True, fast=False)           else:
569           self.__pde_prec.setTolerance(self.getSubProblemTolerance())              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
570           self.__pde_prec.setValue(Y=p)           out=self.__pde_u.getSolution()
571           q=self.__pde_prec.getSolution(verbose=self.show_details)           return  out
572           return q  
573         def norm_Bv(self,Bv):
574        def solve_prec1(self,p):          """
575       #proj=Projector(domain=self.domain, reduce = True, fast=False)          Returns Bv (overwrite).
576           self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
577           self.__pde_prec.setValue(Y=p)          :rtype: equal to the type of p
578           q=self.__pde_prec.getSolution(verbose=self.show_details)          :note: boundary conditions on p should be zero!
579       q0=util.interpolate(q,Function(self.domain))          """
580           print util.inf(q*q0),util.sup(q*q0)          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
581           q-=(1/self.vol)*util.integrate(q0)  
582           print util.inf(q*q0),util.sup(q*q0)       def solve_AinvBt(self,p, tol):
583           return q           """
584             Solves *Av=B^*p* with accuracy `tol`
       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  
           print abs(n_v_old-self.__n_v), n_v, self.getTolerance()  
           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  
585    
586             :param p: a pressure increment
587             :return: the solution of *Av=B^*p*
588             :note: boundary conditions on v should be zero!
589             """
590             self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
591             out=self.__pde_u.getSolution()
592             return  out
593    
594         def solve_prec(self,Bv, tol):
595             """
596             applies preconditioner for for *BA^{-1}B^** to *Bv*
597             with accuracy `self.getSubProblemTolerance()`
598    
599             :param Bv: velocity increment
600             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
601             :note: boundary conditions on p are zero.
602             """
603             self.__pde_prec.setValue(Y=Bv)
604         self.getSolverOptionsPressure().setTolerance(tol)
605         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
606             out=self.__pde_prec.getSolution()
607             return out

Legend:
Removed from v.1562  
changed lines
  Added in v.2881

  ViewVC Help
Powered by ViewVC 1.1.26