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

Contents of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2689 - (show annotations)
Tue Sep 29 05:22:51 2009 UTC (9 years, 10 months ago) by artak
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 23223 byte(s)
getSolverOptionsDiv was returning Pressure options instead of Div equation options
1 ########################################################
2 #
3 # Copyright (c) 2003-2009 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # 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 #
11 ########################################################
12
13 __copyright__="""Copyright (c) 2003-2009 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
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 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37 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):
362 """
363 solves
364
365 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
366 u_{i,i}=0
367
368 u=0 where fixed_u_mask>0
369 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 given 0 is assumed.
372
373 typical usage:
374
375 sp=StokesProblemCartesian(domain)
376 sp.setTolerance()
377 sp.initialize(...)
378 v,p=sp.solve(v0,p0)
379 """
380 def __init__(self,domain,adaptSubTolerance=True, **kwargs):
381 """
382 initialize the Stokes Problem
383
384 :param domain: domain of the problem. The approximation order needs to be two.
385 :type domain: `Domain`
386 :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
387 :type adaptSubTolerance: ``bool``
388 :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
389 """
390 HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
391 self.domain=domain
392 self.vol=util.integrate(1.,Function(self.domain))
393 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
394 self.__pde_u.setSymmetryOn()
395
396 self.__pde_prec=LinearPDE(domain)
397 self.__pde_prec.setReducedOrderOn()
398 self.__pde_prec.setSymmetryOn()
399
400 self.__pde_proj=LinearPDE(domain)
401 self.__pde_proj.setReducedOrderOn()
402 self.__pde_proj.setValue(D=1)
403 self.__pde_proj.setSymmetryOn()
404
405 def getSolverOptionsVelocity(self):
406 """
407 returns the solver options used solve the equation for velocity.
408
409 :rtype: `SolverOptions`
410 """
411 return self.__pde_u.getSolverOptions()
412 def setSolverOptionsVelocity(self, options=None):
413 """
414 set the solver options for solving the equation for velocity.
415
416 :param options: new solver options
417 :type options: `SolverOptions`
418 """
419 self.__pde_u.setSolverOptions(options)
420 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 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 def setSubProblemTolerance(self):
452 """
453 Updates the tolerance for subproblems
454 """
455 if self.adaptSubTolerance():
456 sub_tol=self.getSubProblemTolerance()
457 self.getSolverOptionsDiv().setTolerance(sub_tol)
458 self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
459 self.getSolverOptionsPressure().setTolerance(sub_tol)
460 self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
461 self.getSolverOptionsVelocity().setTolerance(sub_tol)
462 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
463
464
465 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
466 """
467 assigns values to the model parameters
468
469 :param f: external force
470 :type f: `Vector` object in `FunctionSpace` `Function` or similar
471 :param fixed_u_mask: mask of locations with fixed velocity.
472 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
473 :param eta: viscosity
474 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
475 :param surface_stress: normal surface stress
476 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
477 :param stress: initial stress
478 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
479 :note: All values needs to be set.
480 """
481 self.eta=eta
482 A =self.__pde_u.createCoefficient("A")
483 self.__pde_u.setValue(A=Data())
484 for i in range(self.domain.getDim()):
485 for j in range(self.domain.getDim()):
486 A[i,j,j,i] += 1.
487 A[i,j,i,j] += 1.
488 n=self.domain.getNormal()
489 self.__pde_prec.setValue(D=1/self.eta)
490 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
491 self.__f=f
492 self.__surface_stress=surface_stress
493 self.__stress=stress
494
495 def Bv(self,v):
496 """
497 returns inner product of element p and div(v)
498
499 :param v: a residual
500 :return: inner product of element p and div(v)
501 :rtype: ``float``
502 """
503 self.__pde_proj.setValue(Y=-util.div(v))
504 return self.__pde_proj.getSolution()
505
506 def inner_pBv(self,p,Bv):
507 """
508 returns inner product of element p and Bv=-div(v)
509
510 :param p: a pressure increment
511 :param Bv: a residual
512 :return: inner product of element p and Bv=-div(v)
513 :rtype: ``float``
514 """
515 return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
516
517 def inner_p(self,p0,p1):
518 """
519 Returns inner product of p0 and p1
520
521 :param p0: a pressure
522 :param p1: a pressure
523 :return: inner product of p0 and p1
524 :rtype: ``float``
525 """
526 s0=util.interpolate(p0/self.eta,Function(self.domain))
527 s1=util.interpolate(p1/self.eta,Function(self.domain))
528 return util.integrate(s0*s1)
529
530 def norm_v(self,v):
531 """
532 returns the norm of v
533
534 :param v: a velovity
535 :return: norm of v
536 :rtype: non-negative ``float``
537 """
538 return util.sqrt(util.integrate(util.length(util.grad(v))))
539
540 def getV(self, p, v0):
541 """
542 return the value for v for a given p (overwrite)
543
544 :param p: a pressure
545 :param v0: a initial guess for the value v to return.
546 :return: v given as *v= A^{-1} (f-B^*p)*
547 """
548 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
549 if self.__stress.isEmpty():
550 self.__pde_u.setValue(X=p*util.kronecker(self.domain))
551 else:
552 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
553 out=self.__pde_u.getSolution()
554 return out
555
556 def norm_Bv(self,Bv):
557 """
558 Returns Bv (overwrite).
559
560 :rtype: equal to the type of p
561 :note: boundary conditions on p should be zero!
562 """
563 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
564
565 def solve_AinvBt(self,p):
566 """
567 Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`
568
569 :param p: a pressure increment
570 :return: the solution of *Av=B^*p*
571 :note: boundary conditions on v should be zero!
572 """
573 self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
574 out=self.__pde_u.getSolution()
575 return out
576
577 def solve_prec(self,Bv):
578 """
579 applies preconditioner for for *BA^{-1}B^** to *Bv*
580 with accuracy `self.getSubProblemTolerance()`
581
582 :param Bv: velocity increment
583 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
584 :note: boundary conditions on p are zero.
585 """
586 self.__pde_prec.setValue(Y=Bv)
587 return self.__pde_prec.getSolution()

  ViewVC Help
Powered by ViewVC 1.1.26