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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2415 - (hide annotations)
Wed May 13 02:48:39 2009 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 20301 byte(s)
FileWriter added: this class takes care of writing data which are global in MPI to a file. It is recommended to use this class rather then the build in open as it takes care of the case of many processors.
1 ksteube 1809 ########################################################
2 gross 1414 #
3 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
4     # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1414 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1414 #
11 ksteube 1809 ########################################################
12 gross 1414
13 ksteube 1809 __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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1414 """
22     Some models for flow
23    
24     @var __author__: name of author
25     @var __copyright__: copyrights
26     @var __license__: licence agreement
27     @var __url__: url entry point on documentation
28     @var __version__: version
29     @var __date__: date of the version
30     """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36 gross 2100 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38 gross 1414
39 gross 2100 class DarcyFlow(object):
40     """
41 gross 2264 solves the problem
42 gross 1659
43 gross 2100 M{u_i+k_{ij}*p_{,j} = g_i}
44     M{u_{i,i} = f}
45 gross 1659
46 gross 2264 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47 gross 1659
48 gross 2100 @note: The problem is solved in a least squares formulation.
49     """
50 gross 1659
51 gross 2351 def __init__(self, domain, weight=None, useReduced=False):
52 gross 2100 """
53 gross 2208 initializes the Darcy flux problem
54 gross 2100 @param domain: domain of the problem
55     @type domain: L{Domain}
56     """
57     self.domain=domain
58 gross 2351 if weight == None:
59 gross 2354 s=self.domain.getSize()
60     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61 gross 2351 else:
62     self.__l=weight
63 gross 2100 self.__pde_v=LinearPDESystem(domain)
64 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
65     self.__pde_v.setSymmetryOn()
66 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67 gross 2351 # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68 gross 2100 self.__pde_p=LinearSinglePDE(domain)
69     self.__pde_p.setSymmetryOn()
70 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
71 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73 gross 2264 self.setTolerance()
74     self.setAbsoluteTolerance()
75     self.setSubProblemTolerance()
76 gross 1659
77 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
78     """
79 gross 2208 assigns values to model parameters
80 gross 1659
81 gross 2208 @param f: volumetic sources/sinks
82     @type f: scalar value on the domain (e.g. L{Data})
83 gross 2100 @param g: flux sources/sinks
84 gross 2208 @type g: vector values on the domain (e.g. L{Data})
85 gross 2100 @param location_of_fixed_pressure: mask for locations where pressure is fixed
86 gross 2208 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
87     @param location_of_fixed_flux: mask for locations where flux is fixed.
88     @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
89 gross 2264 @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90     C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91 gross 2208 C{v} on the main diagonal is used.
92     @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93 gross 1659
94 gross 2208 @note: the values of parameters which are not set by calling C{setValue} are not altered.
95     @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96     or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
97     is along the M{x_i} axis.
98 gross 2100 """
99 gross 2264 if f !=None:
100 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101     if f.isEmpty():
102     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103     else:
104     if f.getRank()>0: raise ValueError,"illegal rank of f."
105 gross 2267 self.__f=f
106 gross 2264 if g !=None:
107 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108     if g.isEmpty():
109     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110     else:
111     if not g.getShape()==(self.domain.getDim(),):
112     raise ValueError,"illegal shape of g"
113     self.__g=g
114 gross 1659
115 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117 gross 1659
118 gross 2100 if permeability!=None:
119     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120     if perm.getRank()==0:
121     perm=perm*util.kronecker(self.domain.getDim())
122     elif perm.getRank()==1:
123     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125     elif perm.getRank()==2:
126     pass
127     else:
128     raise ValueError,"illegal rank of permeability."
129     self.__permeability=perm
130     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
131 gross 1659
132 gross 2264 def setTolerance(self,rtol=1e-4):
133     """
134     sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135 gross 1659
136 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137    
138     where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
139    
140     @param rtol: relative tolerance for the pressure
141     @type rtol: non-negative C{float}
142 gross 2100 """
143 gross 2264 if rtol<0:
144     raise ValueError,"Relative tolerance needs to be non-negative."
145     self.__rtol=rtol
146     def getTolerance(self):
147 gross 2100 """
148 gross 2264 returns the relative tolerance
149 gross 1659
150 gross 2264 @return: current relative tolerance
151     @rtype: C{float}
152 caltinay 2169 """
153 gross 2264 return self.__rtol
154 gross 1659
155 gross 2264 def setAbsoluteTolerance(self,atol=0.):
156 gross 2208 """
157 gross 2264 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158 gross 1659
159 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160 gross 2208
161 gross 2264 where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
162 gross 2208
163 gross 2264 @param atol: absolute tolerance for the pressure
164     @type atol: non-negative C{float}
165     """
166     if atol<0:
167     raise ValueError,"Absolute tolerance needs to be non-negative."
168     self.__atol=atol
169     def getAbsoluteTolerance(self):
170     """
171     returns the absolute tolerance
172    
173     @return: current absolute tolerance
174     @rtype: C{float}
175     """
176     return self.__atol
177 gross 2208
178 gross 2264 def setSubProblemTolerance(self,rtol=None):
179     """
180     Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181     C{self.getTolerance()**2} is used.
182 gross 2208
183 gross 2264 @param rtol: relative tolerence
184     @type rtol: positive C{float}
185     """
186     if rtol == None:
187     if self.getTolerance()<=0.:
188     raise ValueError,"A positive relative tolerance must be set."
189     self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190     else:
191     if rtol<=0:
192     raise ValueError,"sub-problem tolerance must be positive."
193     self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194 gross 2208
195 gross 2264 def getSubProblemTolerance(self):
196     """
197     Returns the subproblem reduction factor.
198 gross 2208
199 gross 2264 @return: subproblem reduction factor
200     @rtype: C{float}
201     """
202     return self.__sub_tol
203 gross 2208
204 gross 2264 def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205     """
206     solves the problem.
207 gross 1659
208 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
209 gross 1659
210 gross 2208 @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
211     @type u0: vector value on the domain (e.g. L{Data}).
212     @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
213     @type p0: scalar value on the domain (e.g. L{Data}).
214     @param verbose: if set some information on iteration progress are printed
215     @type verbose: C{bool}
216     @param show_details: if set information on the subiteration process are printed.
217     @type show_details: C{bool}
218     @return: flux and pressure
219     @rtype: C{tuple} of L{Data}.
220 gross 2264
221 gross 2208 @note: The problem is solved as a least squares form
222 gross 2100
223 gross 2208 M{(I+D^*D)u+Qp=D^*f+g}
224     M{Q^*u+Q^*Qp=Q^*g}
225 gross 2100
226 gross 2264 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227     We eliminate the flux form the problem by setting
228 caltinay 2169
229 gross 2208 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
230 caltinay 2169
231 gross 2208 form the first equation. Inserted into the second equation we get
232 caltinay 2169
233 gross 2208 M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0 on location_of_fixed_pressure
234 gross 2264
235     which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
236 gross 2208 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237     """
238 gross 2386 self.verbose=verbose
239 gross 2208 self.show_details= show_details and self.verbose
240 gross 2264 rtol=self.getTolerance()
241     atol=self.getAbsoluteTolerance()
242     if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243 caltinay 2169
244 gross 2264 num_corrections=0
245     converged=False
246     p=p0
247     norm_r=None
248     while not converged:
249     v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250     Qp=self.__Q(p)
251     norm_v=self.__L2(v)
252     norm_Qp=self.__L2(Qp)
253     if norm_v == 0.:
254     if norm_Qp == 0.:
255     return v,p
256     else:
257     fac=norm_Qp
258     else:
259     if norm_Qp == 0.:
260     fac=norm_v
261     else:
262     fac=2./(1./norm_v+1./norm_Qp)
263     ATOL=(atol+rtol*fac)
264     if self.verbose:
265     print "DarcyFlux: L2 norm of v = %e."%norm_v
266     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268     if norm_r == None or norm_r>ATOL:
269     if num_corrections>max_num_corrections:
270     raise ValueError,"maximum number of correction steps reached."
271 gross 2354 p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272 gross 2264 num_corrections+=1
273     else:
274     converged=True
275     return v,p
276     def __L2(self,v):
277     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278    
279     def __Q(self,p):
280     return util.tensor_mult(self.__permeability,util.grad(p))
281    
282     def __Aprod(self,dp):
283     self.__pde_v.setTolerance(self.getSubProblemTolerance())
284 gross 2208 if self.show_details: print "DarcyFlux: Applying operator"
285 gross 2264 Qdp=self.__Q(dp)
286     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287 gross 2349 du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
289 gross 2264 return Qdp+du
290     def __inner_GMRES(self,r,s):
291     return util.integrate(util.inner(r,s))
292    
293 gross 2100 def __inner_PCG(self,p,r):
294 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
295 gross 2100
296     def __Msolve_PCG(self,r):
297 gross 2264 self.__pde_p.setTolerance(self.getSubProblemTolerance())
298 gross 2208 if self.show_details: print "DarcyFlux: Applying preconditioner"
299 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
301 gross 2349 return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302 gross 2100
303 gross 2264 def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304     """
305     returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
306     on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
307     Note that C{g} and C{f} are used, see L{setValue}.
308    
309     @param p: pressure.
310     @type p: scalar value on the domain (e.g. L{Data}).
311     @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
312     @type fixed_flux: vector values on the domain (e.g. L{Data}).
313     @param tol: relative tolerance to be used.
314     @type tol: positive C{float}.
315     @return: flux
316     @rtype: L{Data}
317     @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
318     for the permeability M{k_{ij}}
319     """
320     self.__pde_v.setTolerance(self.getSubProblemTolerance())
321     g=self.__g
322     f=self.__f
323 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
324 gross 2264 if p == None:
325     self.__pde_v.setValue(Y=g)
326     else:
327     self.__pde_v.setValue(Y=g-self.__Q(p))
328 gross 2349 return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329 gross 2264
330 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
331 gross 2251 """
332 gross 2264 solves
333 gross 1414
334 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
335     u_{i,i}=0
336 gross 1414
337 gross 2208 u=0 where fixed_u_mask>0
338     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
339 gross 1414
340 gross 2264 if surface_stress is not given 0 is assumed.
341 gross 1414
342 gross 2251 typical usage:
343 gross 1414
344 gross 2208 sp=StokesProblemCartesian(domain)
345     sp.setTolerance()
346     sp.initialize(...)
347     v,p=sp.solve(v0,p0)
348 gross 2251 """
349     def __init__(self,domain,**kwargs):
350 gross 2100 """
351 gross 2208 initialize the Stokes Problem
352 gross 2100
353 gross 2208 @param domain: domain of the problem. The approximation order needs to be two.
354 gross 2100 @type domain: L{Domain}
355 gross 2208 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
356 gross 2100 """
357 gross 1414 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
358     self.domain=domain
359     self.vol=util.integrate(1.,Function(self.domain))
360     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
361     self.__pde_u.setSymmetryOn()
362 gross 2100 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363 gross 2351 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364 gross 2264
365 gross 1414 self.__pde_prec=LinearPDE(domain)
366     self.__pde_prec.setReducedOrderOn()
367 gross 2156 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368 gross 1414 self.__pde_prec.setSymmetryOn()
369    
370 gross 2415 self.__pde_proj=LinearPDE(domain)
371     self.__pde_proj.setReducedOrderOn()
372     self.__pde_proj.setValue(D=1)
373     self.__pde_proj.setSymmetryOn()
374    
375    
376 gross 2251 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
377 gross 2208 """
378     assigns values to the model parameters
379 gross 2100
380 gross 2208 @param f: external force
381     @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
382     @param fixed_u_mask: mask of locations with fixed velocity.
383     @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
384     @param eta: viscosity
385     @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
386     @param surface_stress: normal surface stress
387     @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
388     @param stress: initial stress
389     @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
390     @note: All values needs to be set.
391 gross 2264
392 gross 2208 """
393     self.eta=eta
394     A =self.__pde_u.createCoefficient("A")
395     self.__pde_u.setValue(A=Data())
396     for i in range(self.domain.getDim()):
397     for j in range(self.domain.getDim()):
398 gross 2264 A[i,j,j,i] += 1.
399 gross 2208 A[i,j,i,j] += 1.
400 gross 2264 self.__pde_prec.setValue(D=1/self.eta)
401 gross 2251 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)
402     self.__f=f
403     self.__surface_stress=surface_stress
404 gross 2208 self.__stress=stress
405 gross 1414
406 gross 2251 def inner_pBv(self,p,v):
407     """
408     returns inner product of element p and div(v)
409 gross 1414
410 gross 2251 @param p: a pressure increment
411     @param v: a residual
412     @return: inner product of element p and div(v)
413     @rtype: C{float}
414 gross 2100 """
415 gross 2251 return util.integrate(-p*util.div(v))
416 gross 2208
417 gross 2251 def inner_p(self,p0,p1):
418 gross 2100 """
419 gross 2251 Returns inner product of p0 and p1
420 gross 1414
421 gross 2251 @param p0: a pressure
422     @param p1: a pressure
423     @return: inner product of p0 and p1
424 gross 2208 @rtype: C{float}
425 gross 2100 """
426     s0=util.interpolate(p0/self.eta,Function(self.domain))
427     s1=util.interpolate(p1/self.eta,Function(self.domain))
428     return util.integrate(s0*s1)
429 artak 1550
430 gross 2251 def norm_v(self,v):
431 gross 2100 """
432 gross 2251 returns the norm of v
433 gross 2208
434 gross 2251 @param v: a velovity
435     @return: norm of v
436     @rtype: non-negative C{float}
437 gross 2100 """
438 gross 2251 return util.sqrt(util.integrate(util.length(util.grad(v))))
439 gross 2100
440 gross 2251 def getV(self, p, v0):
441 gross 1414 """
442 gross 2251 return the value for v for a given p (overwrite)
443    
444     @param p: a pressure
445 gross 2264 @param v0: a initial guess for the value v to return.
446 gross 2251 @return: v given as M{v= A^{-1} (f-B^*p)}
447 gross 1414 """
448     self.__pde_u.setTolerance(self.getSubProblemTolerance())
449 gross 2251 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
450 gross 2100 if self.__stress.isEmpty():
451 gross 2251 self.__pde_u.setValue(X=p*util.kronecker(self.domain))
452 gross 2100 else:
453 gross 2251 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
454 gross 2100 out=self.__pde_u.getSolution(verbose=self.show_details)
455 gross 2208 return out
456 gross 1414
457 gross 2251
458     raise NotImplementedError,"no v calculation implemented."
459    
460 gross 2264
461 gross 2251 def norm_Bv(self,v):
462     """
463     Returns Bv (overwrite).
464    
465     @rtype: equal to the type of p
466     @note: boundary conditions on p should be zero!
467     """
468 gross 2415 self.__pde_proj.setValue(Y=util.div(v))
469     self.__pde_prec.setTolerance(self.getSubProblemTolerance())
470     return util.sqrt(util.integrate(util.interpolate(self.__pde_proj.getSolution(),Function(self.domain))**2))
471 gross 2251
472     def solve_AinvBt(self,p):
473     """
474     Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
475    
476     @param p: a pressure increment
477 gross 2264 @return: the solution of M{Av=B^*p}
478 gross 2251 @note: boundary conditions on v should be zero!
479     """
480     self.__pde_u.setTolerance(self.getSubProblemTolerance())
481     self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
482     out=self.__pde_u.getSolution(verbose=self.show_details)
483     return out
484    
485     def solve_precB(self,v):
486     """
487     applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
488     with accuracy L{self.getSubProblemTolerance()} (overwrite).
489    
490     @param v: velocity increment
491     @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
492     @note: boundary conditions on p are zero.
493     """
494     self.__pde_prec.setValue(Y=-util.div(v))
495 gross 1414 self.__pde_prec.setTolerance(self.getSubProblemTolerance())
496 gross 2251 return self.__pde_prec.getSolution(verbose=self.show_details)

  ViewVC Help
Powered by ViewVC 1.1.26