/[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 1467 by gross, Wed Apr 2 08:10:37 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 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__="http://www.uq.edu.au/esscc/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  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
38    
39    class DarcyFlow(object):
40        """
41        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,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            self.__pde_v=LinearPDESystem(domain)
59            if useReduced: self.__pde_v.setReducedOrderOn()
60            self.__pde_v.setSymmetryOn()
61            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
62            self.__pde_p=LinearSinglePDE(domain)
63            self.__pde_p.setSymmetryOn()
64            if useReduced: self.__pde_p.setReducedOrderOn()
65            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67            self.__ATOL= None
68    
69        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70            """
71            assigns values to model parameters
72    
73            @param f: volumetic sources/sinks
74            @type f: scalar value on the domain (e.g. L{Data})
75            @param g: flux sources/sinks
76            @type g: vector values on the domain (e.g. L{Data})
77            @param location_of_fixed_pressure: mask for locations where pressure is fixed
78            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
79            @param location_of_fixed_flux:  mask for locations where flux is fixed.
80            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
81            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
82                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
83                                 C{v} on the main diagonal is used.
84            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
85    
86            @note: the values of parameters which are not set by calling C{setValue} are not altered.
87            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
88                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
89                   is along the M{x_i} axis.
90            """
91            if f !=None:
92               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
93               if f.isEmpty():
94                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
95               else:
96                   if f.getRank()>0: raise ValueError,"illegal rank of f."
97               self.f=f
98            if g !=None:  
99               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
100               if g.isEmpty():
101                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
102               else:
103                 if not g.getShape()==(self.domain.getDim(),):
104                   raise ValueError,"illegal shape of g"
105               self.__g=g
106    
107            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
108            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
109    
110            if permeability!=None:
111               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
112               if perm.getRank()==0:
113                   perm=perm*util.kronecker(self.domain.getDim())
114               elif perm.getRank()==1:
115                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
116                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
117               elif perm.getRank()==2:
118                  pass
119               else:
120                  raise ValueError,"illegal rank of permeability."
121               self.__permeability=perm
122               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
123    
124    
125        def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
126            """
127            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
128            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
129            Note that C{g} and C{f} are used, see L{setValue}.
130            
131            @param p: pressure.
132            @type p: scalar value on the domain (e.g. L{Data}).
133            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134            @type fixed_flux: vector values on the domain (e.g. L{Data}).
135            @param tol: relative tolerance to be used.
136            @type tol: positive C{float}.
137            @return: flux
138            @rtype: L{Data}
139            @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}}
140                   for the permeability M{k_{ij}}
141            """
142            self.__pde_v.setTolerance(tol)
143            g=self.__g
144            f=self.__f
145            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
146            if p == None:
147               self.__pde_v.setValue(Y=g)
148            else:
149               self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
150            return self.__pde_v.getSolution(verbose=show_details)
151    
152        def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
153            """
154            returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
155            on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).
156            Note that C{g} is used, see L{setValue}.
157            
158            @param v: flux.
159            @type v: vector-valued on the domain (e.g. L{Data}).
160            @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
161            @type fixed_pressure: vector values on the domain (e.g. L{Data}).
162            @param tol: relative tolerance to be used.
163            @type tol: positive C{float}.
164            @return: pressure
165            @rtype: L{Data}
166            @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}
167                   for the permeability M{k_{ij}}
168            """
169            self.__pde_v.setTolerance(tol)
170            g=self.__g
171            self.__pde_p.setValue(r=fixed_pressure)
172            if v == None:
173               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
174            else:
175               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
176            return self.__pde_p.getSolution(verbose=show_details)
177    
178        def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
179            """
180            set the tolerance C{ATOL} used to terminate the solution process. It is used
181    
182            M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}
183    
184            where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.
185    
186            The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual
187    
188            M{r=Q^*(g-Qp-v)}
189    
190            the condition
191    
192            M{<(Q^*Q)^{-1} r,r> <= ATOL}
193    
194            holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
195    
196            @param atol: absolute tolerance for the pressure
197            @type atol: non-negative C{float}
198            @param rtol: relative tolerance for the pressure
199            @type rtol: non-negative C{float}
200            @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the
201            L{getPressure} method or use  the value from a previous time step.
202            @type p_ref: scalar value on the domain (e.g. L{Data}).
203            @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the
204            L{getFlux} method or use  the value from a previous time step.
205            @type v_ref: vector-valued on the domain (e.g. L{Data}).
206            @return: used absolute tolerance.
207            @rtype: positive C{float}
208            """
209            g=self.__g
210            if not v_ref == None:
211               f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212            else:
213               f1=util.integrate(util.length(util.interpolate(g))**2)
214            if not p_ref == None:
215               f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
216            else:
217               f2=0
218            self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219            if self.__ATOL<=0:
220               raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221            return self.__ATOL
222            
223        def getTolerance(self):
224            """
225            returns the current tolerance.
226      
227            @return: used absolute tolerance.
228            @rtype: positive C{float}
229            """
230            if self.__ATOL==None:
231               raise ValueError,"no tolerance is defined."
232            return self.__ATOL
233    
234        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235             """
236             solves the problem.
237    
238             The iteration is terminated if the residual norm is less then self.getTolerance().
239    
240             @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.
241             @type u0: vector value on the domain (e.g. L{Data}).
242             @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.
243             @type p0: scalar value on the domain (e.g. L{Data}).
244             @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245             @type sub_rtol: positive-negative C{float}
246             @param verbose: if set some information on iteration progress are printed
247             @type verbose: C{bool}
248             @param show_details:  if set information on the subiteration process are printed.
249             @type show_details: C{bool}
250             @return: flux and pressure
251             @rtype: C{tuple} of L{Data}.
252    
253             @note: The problem is solved as a least squares form
254    
255             M{(I+D^*D)u+Qp=D^*f+g}
256             M{Q^*u+Q^*Qp=Q^*g}
257    
258             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
259             We eliminate the flux form the problem by setting
260    
261             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
262    
263             form the first equation. Inserted into the second equation we get
264    
265             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
266            
267             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
268             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
269             """
270             self.verbose=verbose
271             self.show_details= show_details and self.verbose
272             self.__pde_v.setTolerance(sub_rtol)
273             self.__pde_p.setTolerance(sub_rtol)
274             ATOL=self.getTolerance()
275             if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276             #########################################################################################################################
277             #
278             #   we solve:
279             #  
280             #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
281             #
282             #   residual is
283             #
284             #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)
285             #
286             #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
287             #
288             #    we use (g - Qp, v) to represent the residual. not that
289             #
290             #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291             #
292             #   while the initial residual is
293             #
294             #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295             #  
296             d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297             self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
298             v00=self.__pde_v.getSolution(verbose=show_details)
299             if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300             self.__pde_v.setValue(r=Data())
301             # start CG
302             r=ArithmeticTuple(d0, v00)
303             p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
304             return r[1],p
305    
306        def __Aprod_PCG(self,dp):
307              if self.show_details: print "DarcyFlux: Applying operator"
308              #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
309              mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
310              self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311              du=self.__pde_v.getSolution(verbose=self.show_details)
312              return ArithmeticTuple(mQdp,du)
313    
314        def __inner_PCG(self,p,r):
315             a=util.tensor_mult(self.__permeability,util.grad(p))
316             f0=util.integrate(util.inner(a,r[0]))
317             f1=util.integrate(util.inner(a,r[1]))
318             # print "__inner_PCG:",f0,f1,"->",f0-f1
319             return f0-f1
320    
321        def __Msolve_PCG(self,r):
322              if self.show_details: print "DarcyFlux: Applying preconditioner"
323              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
324              return self.__pde_p.getSolution(verbose=self.show_details)
325    
326  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
327        """        """
328        solves        solves
329    
330            -(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}
331                  u_{i,i}=0                  u_{i,i}=0
332    
333            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
334            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
335    
336        if surface_stress is not give 0 is assumed.        if surface_stress is not given 0 is assumed.
337    
338        typical usage:        typical usage:
339    
# Line 58  class StokesProblemCartesian(Homogeneous Line 343  class StokesProblemCartesian(Homogeneous
343              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
344        """        """
345        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
346             """
347             initialize the Stokes Problem
348    
349             @param domain: domain of the problem. The approximation order needs to be two.
350             @type domain: L{Domain}
351             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
352             """
353           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
354           self.domain=domain           self.domain=domain
355           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
356           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())
357           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
358           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
359             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
360                            
361           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
362           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
363             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
365    
366           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 74  class StokesProblemCartesian(Homogeneous Line 368  class StokesProblemCartesian(Homogeneous
368           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
369           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
370    
371        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
372            """
373            assigns values to the model parameters
374    
375            @param f: external force
376            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
377            @param fixed_u_mask: mask of locations with fixed velocity.
378            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
379            @param eta: viscosity
380            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
381            @param surface_stress: normal surface stress
382            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
383            @param stress: initial stress
384        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
385            @note: All values needs to be set.
386    
387            """
388          self.eta=eta          self.eta=eta
389          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
390      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
391          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
392          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
393              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
394              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
395      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
396          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,Y=f,y=surface_stress)
397            self.__stress=stress
398    
399          def B(self,v):
400            """
401            returns div(v)
402            @rtype: equal to the type of p
403    
404            @note: boundary conditions on p should be zero!
405            """
406            if self.show_details: print "apply divergence:"
407            self.__pde_proj.setValue(Y=-util.div(v))
408            self.__pde_proj.setTolerance(self.getSubProblemTolerance())
409            return self.__pde_proj.getSolution(verbose=self.show_details)
410    
411        def B(self,arg):        def inner_pBv(self,p,Bv):
412           d=util.div(arg)           """
413           self.__pde_proj.setValue(Y=d)           returns inner product of element p and Bv  (overwrite)
414           self.__pde_proj.setTolerance(self.getSubProblemTolerance())          
415           return self.__pde_proj.getSolution(verbose=self.show_details)           @type p: equal to the type of p
416             @type Bv: equal to the type of result of operator B
417        def inner(self,p0,p1):           @rtype: C{float}
418           s0=util.interpolate(p0,Function(self.domain))  
419           s1=util.interpolate(p1,Function(self.domain))           @rtype: equal to the type of p
420             """
421             s0=util.interpolate(p,Function(self.domain))
422             s1=util.interpolate(Bv,Function(self.domain))
423           return util.integrate(s0*s1)           return util.integrate(s0*s1)
424    
425        def getStress(self,u):        def inner_p(self,p0,p1):
426           mg=util.grad(u)           """
427           return 2.*self.eta*util.symmetric(mg)           returns inner product of element p0 and p1  (overwrite)
428            
429             @type p0: equal to the type of p
430             @type p1: equal to the type of p
431             @rtype: C{float}
432    
433             @rtype: equal to the type of p
434             """
435             s0=util.interpolate(p0/self.eta,Function(self.domain))
436             s1=util.interpolate(p1/self.eta,Function(self.domain))
437             return util.integrate(s0*s1)
438    
439          def inner_v(self,v0,v1):
440             """
441             returns inner product of two element v0 and v1  (overwrite)
442            
443             @type v0: equal to the type of v
444             @type v1: equal to the type of v
445             @rtype: C{float}
446    
447             @rtype: equal to the type of v
448             """
449         gv0=util.grad(v0)
450         gv1=util.grad(v1)
451             return util.integrate(util.inner(gv0,gv1))
452    
453        def solve_A(self,u,p):        def solve_A(self,u,p):
454           """           """
455           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
456           """           """
457             if self.show_details: print "solve for velocity:"
458           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
459           self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain))           if self.__stress.isEmpty():
460           return  self.__pde_u.getSolution(verbose=self.show_details)              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
461             else:
462                self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
463             out=self.__pde_u.getSolution(verbose=self.show_details)
464             return  out
465    
466        def solve_prec(self,p):        def solve_prec(self,p):
467             if self.show_details: print "apply preconditioner:"
468           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
469           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
470           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
471           return q           return q
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.Lsup(v)  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)  
           self.iter+=1  
           if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium_GMRES(self,norm_r,norm_b):  
           if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())  
           self.iter+=1  
           if norm_r <= norm_b*self.getTolerance():  
               if self.verbose: print "GMRES terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  

Legend:
Removed from v.1467  
changed lines
  Added in v.2208

  ViewVC Help
Powered by ViewVC 1.1.26