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

revision 1562 by gross, Wed May 21 13:04:40 2008 UTC revision 2288 by gross, Tue Feb 24 06:11:48 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
#
#######################################################
10  #  #
11    ########################################################
12
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
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"
http://www.access.edu.au
__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 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.__l=util.longestEdge(self.domain)**2
59            self.__pde_v=LinearPDESystem(domain)
60            if useReduced: self.__pde_v.setReducedOrderOn()
61            self.__pde_v.setSymmetryOn()
62            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
63            self.__pde_p=LinearSinglePDE(domain)
64            self.__pde_p.setSymmetryOn()
65            if useReduced: self.__pde_p.setReducedOrderOn()
66            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
67            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
68            self.setTolerance()
69            self.setAbsoluteTolerance()
70            self.setSubProblemTolerance()
71
72        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
73            """
74            assigns values to model parameters
75
76            @param f: volumetic sources/sinks
77            @type f: scalar value on the domain (e.g. L{Data})
78            @param g: flux sources/sinks
79            @type g: vector values on the domain (e.g. L{Data})
80            @param location_of_fixed_pressure: mask for locations where pressure is fixed
81            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
82            @param location_of_fixed_flux:  mask for locations where flux is fixed.
83            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
84            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
85                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
86                                 C{v} on the main diagonal is used.
87            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
88
89            @note: the values of parameters which are not set by calling C{setValue} are not altered.
90            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
91                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
92                   is along the M{x_i} axis.
93            """
94            if f !=None:
95               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
96               if f.isEmpty():
97                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
98               else:
99                   if f.getRank()>0: raise ValueError,"illegal rank of f."
100               self.__f=f
101            if g !=None:
102               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
103               if g.isEmpty():
104                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
105               else:
106                 if not g.getShape()==(self.domain.getDim(),):
107                   raise ValueError,"illegal shape of g"
108               self.__g=g
109
110            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
111            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
112
113            if permeability!=None:
114               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
115               if perm.getRank()==0:
116                   perm=perm*util.kronecker(self.domain.getDim())
117               elif perm.getRank()==1:
118                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
119                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
120               elif perm.getRank()==2:
121                  pass
122               else:
123                  raise ValueError,"illegal rank of permeability."
124               self.__permeability=perm
125               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
126
127        def setTolerance(self,rtol=1e-4):
128            """
129            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
130
131            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
132
133            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}}.
134
135            @param rtol: relative tolerance for the pressure
136            @type rtol: non-negative C{float}
137            """
138            if rtol<0:
139                raise ValueError,"Relative tolerance needs to be non-negative."
140            self.__rtol=rtol
141        def getTolerance(self):
142            """
143            returns the relative tolerance
144
145            @return: current relative tolerance
146            @rtype: C{float}
147            """
148            return self.__rtol
149
150        def setAbsoluteTolerance(self,atol=0.):
151            """
152            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
153
154            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
155
156            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}}.
157
158            @param atol: absolute tolerance for the pressure
159            @type atol: non-negative C{float}
160            """
161            if atol<0:
162                raise ValueError,"Absolute tolerance needs to be non-negative."
163            self.__atol=atol
164        def getAbsoluteTolerance(self):
165           """
166           returns the absolute tolerance
167
168           @return: current absolute tolerance
169           @rtype: C{float}
170           """
171           return self.__atol
172
173        def setSubProblemTolerance(self,rtol=None):
174             """
175             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
176             C{self.getTolerance()**2} is used.
177
178             @param rtol: relative tolerence
179             @type rtol: positive C{float}
180             """
181             if rtol == None:
182                  if self.getTolerance()<=0.:
183                      raise ValueError,"A positive relative tolerance must be set."
184                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
185             else:
186                 if rtol<=0:
187                     raise ValueError,"sub-problem tolerance must be positive."
188                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
189
190        def getSubProblemTolerance(self):
191             """
192             Returns the subproblem reduction factor.
193
194             @return: subproblem reduction factor
195             @rtype: C{float}
196             """
197             return self.__sub_tol
198
199        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
200             """
201             solves the problem.
202
203             The iteration is terminated if the residual norm is less then self.getTolerance().
204
205             @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.
206             @type u0: vector value on the domain (e.g. L{Data}).
207             @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.
208             @type p0: scalar value on the domain (e.g. L{Data}).
209             @param verbose: if set some information on iteration progress are printed
210             @type verbose: C{bool}
211             @param show_details:  if set information on the subiteration process are printed.
212             @type show_details: C{bool}
213             @return: flux and pressure
214             @rtype: C{tuple} of L{Data}.
215
216             @note: The problem is solved as a least squares form
217
218             M{(I+D^*D)u+Qp=D^*f+g}
219             M{Q^*u+Q^*Qp=Q^*g}
220
221             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
222             We eliminate the flux form the problem by setting
223
224             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
225
226             form the first equation. Inserted into the second equation we get
227
228             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
229
230             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
231             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
232             """
233             self.verbose=verbose or True
234             self.show_details= show_details and self.verbose
235             rtol=self.getTolerance()
236             atol=self.getAbsoluteTolerance()
237             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
238
239             num_corrections=0
240             converged=False
241             p=p0
242             norm_r=None
243             while not converged:
244                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
245                   Qp=self.__Q(p)
246                   norm_v=self.__L2(v)
247                   norm_Qp=self.__L2(Qp)
248                   if norm_v == 0.:
249                      if norm_Qp == 0.:
250                         return v,p
251                      else:
252                        fac=norm_Qp
253                   else:
254                      if norm_Qp == 0.:
255                        fac=norm_v
256                      else:
257                        fac=2./(1./norm_v+1./norm_Qp)
258                   ATOL=(atol+rtol*fac)
259                   if self.verbose:
260                        print "DarcyFlux: L2 norm of v = %e."%norm_v
261                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
262                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
263                   if norm_r == None or norm_r>ATOL:
264                       if num_corrections>max_num_corrections:
265                             raise ValueError,"maximum number of correction steps reached."
266                       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.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
267                       num_corrections+=1
268                   else:
269                       converged=True
270             return v,p
271    #
272    #
273    #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp
274    #               #===========================================================================
275    #               norm_r_hat=self.__L2(r_hat)
276    #               norm_v=self.__L2(v)
277    #               norm_g=self.__L2(g)
278    #               norm_gv=self.__L2(g-v)
279    #               norm_Qp=self.__L2(Qp)
280    #               norm_gQp=self.__L2(g-Qp)
281    #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))
282    #               fac=min(norm_v,norm_Qp,norm_gv)
283    #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
284    #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat
285    #               if r!=None:
286    #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat
287    #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()
288    #                   self.setSubProblemTolerance(sub_tol)
289    #                   print "subtol_new=",self.getSubProblemTolerance()
290    #               print "norm_v = ",norm_v
291    #               print "norm_gv = ",norm_gv
292    #               print "norm_Qp = ",norm_Qp
293    #               print "norm_gQp = ",norm_gQp
294    #               print "norm_g = ",norm_g
295    #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)
296    #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)
297    #               if fac == 0:
298    #                   if self.verbose: print "DarcyFlux: trivial case!"
299    #                   return v,p
300    #               #===============================================================================
301    #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))
302    #               # norm_Qp=self.__L2(Qp)
303    #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))
304    #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp
305    #
306    #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))
307    #               if self.verbose:
308    #                   print "DarcyFlux: residual = %e"%norm_r_hat
309    #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
310    #               if norm_r_hat <= ATOL:
311    #                   print "DarcyFlux: iteration finalized."
312    #                   converged=True
313    #               else:
314    #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)
315    #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)
316    #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
317    #               print "norm_r =",norm_r
318    #         return v,p
319        def __L2(self,v):
320             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
321
322        def __Q(self,p):
324
325        def __Aprod(self,dp):
326              self.__pde_v.setTolerance(self.getSubProblemTolerance())
327              if self.show_details: print "DarcyFlux: Applying operator"
328              Qdp=self.__Q(dp)
329              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
330              du=self.__pde_v.getSolution(verbose=self.show_details)
331              return Qdp+du
332        def __inner_GMRES(self,r,s):
333             return util.integrate(util.inner(r,s))
334
335        def __inner_PCG(self,p,r):
336             return util.integrate(util.inner(self.__Q(p), r))
337
338        def __Msolve_PCG(self,r):
339              self.__pde_p.setTolerance(self.getSubProblemTolerance())
340              if self.show_details: print "DarcyFlux: Applying preconditioner"
341              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
342              return self.__pde_p.getSolution(verbose=self.show_details)
343
344        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
345            """
346            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
347            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
348            Note that C{g} and C{f} are used, see L{setValue}.
349
350            @param p: pressure.
351            @type p: scalar value on the domain (e.g. L{Data}).
352            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
353            @type fixed_flux: vector values on the domain (e.g. L{Data}).
354            @param tol: relative tolerance to be used.
355            @type tol: positive C{float}.
356            @return: flux
357            @rtype: L{Data}
358            @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}}
359                   for the permeability M{k_{ij}}
360            """
361            self.__pde_v.setTolerance(self.getSubProblemTolerance())
362            g=self.__g
363            f=self.__f
364            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
365            if p == None:
366               self.__pde_v.setValue(Y=g)
367            else:
368               self.__pde_v.setValue(Y=g-self.__Q(p))
369            return self.__pde_v.getSolution(verbose=show_details)
370
372        """       """
373        solves       solves
374
375            -(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}
376                  u_{i,i}=0                  u_{i,i}=0
377
379            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
380
381        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
382
383        typical usage:       typical usage:
384
385              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
386              sp.setTolerance()              sp.setTolerance()
387              sp.initialize(...)              sp.initialize(...)
388              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
389        """       """
390        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
391             """
392             initialize the Stokes Problem
393
394             @param domain: domain of the problem. The approximation order needs to be two.
395             @type domain: L{Domain}
396             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
397             """
399           self.domain=domain           self.domain=domain
400           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
401           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())
402           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
403           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
404                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
405
406           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
407           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
408             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
409           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
410
412           self.__pde_proj.setReducedOrderOn()          """
413           self.__pde_proj.setSymmetryOn()          assigns values to the model parameters
414           self.__pde_proj.setValue(D=1.)
415            @param f: external force
416            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
418            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
419            @param eta: viscosity
420            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
421            @param surface_stress: normal surface stress
422            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
423            @param stress: initial stress
424        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
425            @note: All values needs to be set.
426
428          self.eta=eta          self.eta=eta
429          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
430      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
431          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
432          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
433              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
434              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
435      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
437            self.__f=f
438            self.__surface_stress=surface_stress
439            self.__stress=stress
440
441         def inner_pBv(self,p,v):
442             """
443             returns inner product of element p and div(v)
444
445             @param p: a pressure increment
446             @param v: a residual
447             @return: inner product of element p and div(v)
448             @rtype: C{float}
449             """
450             return util.integrate(-p*util.div(v))
451
452         def inner_p(self,p0,p1):
453             """
454             Returns inner product of p0 and p1
455
456        def B(self,arg):           @param p0: a pressure
457           d=util.div(arg)           @param p1: a pressure
458           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
459           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
460           return self.__pde_proj.getSolution(verbose=self.show_details)           """
461             s0=util.interpolate(p0/self.eta,Function(self.domain))
462        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))
463           return util.integrate(s0*s1)           return util.integrate(s0*s1)
464
465        def inner_a(self,a0,a1):       def norm_v(self,v):
466           p0=util.interpolate(a0[1],Function(self.domain))           """
467           p1=util.interpolate(a1[1],Function(self.domain))           returns the norm of v
alfa=(1/self.vol)*util.integrate(p0)
beta=(1/self.vol)*util.integrate(p1)
#print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))
return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))

def getStress(self,u):
return 2.*self.eta*util.symmetric(mg)
468
469        def solve_A(self,u,p):           @param v: a velovity
470             @return: norm of v
471             @rtype: non-negative C{float}
472           """           """
474
475         def getV(self, p, v0):
476             """
477             return the value for v for a given p (overwrite)
478
479             @param p: a pressure
480             @param v0: a initial guess for the value v to return.
481             @return: v given as M{v= A^{-1} (f-B^*p)}
482           """           """
483           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
484           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)
485           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
486                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
487             else:
488                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
489             out=self.__pde_u.getSolution(verbose=self.show_details)
490             return  out
491
492
493        def solve_prec(self,p):           raise NotImplementedError,"no v calculation implemented."
#proj=Projector(domain=self.domain, reduce = True, fast=False)
self.__pde_prec.setTolerance(self.getSubProblemTolerance())
self.__pde_prec.setValue(Y=p)
q=self.__pde_prec.getSolution(verbose=self.show_details)
return q
494
def solve_prec1(self,p):
#proj=Projector(domain=self.domain, reduce = True, fast=False)
self.__pde_prec.setTolerance(self.getSubProblemTolerance())
self.__pde_prec.setValue(Y=p)
q=self.__pde_prec.getSolution(verbose=self.show_details)
q0=util.interpolate(q,Function(self.domain))
print util.inf(q*q0),util.sup(q*q0)
q-=(1/self.vol)*util.integrate(q0)
print util.inf(q*q0),util.sup(q*q0)
return q

def stoppingcriterium(self,Bv,v,p):
n_r=util.sqrt(self.inner(Bv,Bv))
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
495
496         def norm_Bv(self,v):
497            """
498            Returns Bv (overwrite).
499
500            @rtype: equal to the type of p
501            @note: boundary conditions on p should be zero!
502            """
503            return util.sqrt(util.integrate(util.div(v)**2))
504
505         def solve_AinvBt(self,p):
506             """
507             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
508
509             @param p: a pressure increment
510             @return: the solution of M{Av=B^*p}
511             @note: boundary conditions on v should be zero!
512             """
513             self.__pde_u.setTolerance(self.getSubProblemTolerance())
514             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
515             out=self.__pde_u.getSolution(verbose=self.show_details)
516             return  out
517
518         def solve_precB(self,v):
519             """
520             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
521             with accuracy L{self.getSubProblemTolerance()} (overwrite).
522
523             @param v: velocity increment
524             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
525             @note: boundary conditions on p are zero.
526             """
527             self.__pde_prec.setValue(Y=-util.div(v))
528             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
529             return self.__pde_prec.getSolution(verbose=self.show_details)

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