# Contents of /trunk/escript/py_src/flows.py

Revision 2486 - (show annotations)
Tue Jun 23 03:38:54 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 23295 byte(s)
```a bit more info are printed for the DarcyFlux now.
```
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2008 by University of Queensland 4 # Earth Systems Science Computational Center (ESSCC) 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 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 19 __url__= 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 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, weight=None, useReduced=False, adaptSubTolerance=True): 52 """ 53 initializes the Darcy flux problem 54 @param domain: domain of the problem 55 @type domain: L{Domain} 56 @param useReduced: uses reduced oreder on flux and pressure 57 @type useReduced: C{bool} 58 @param adaptSubTolerance: switches on automatic subtolerance selection 59 @type adaptSubTolerance: C{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 M{(I+D^*D)u=F} 87 88 @return: L{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 M{(I+D^*D)u=F} 96 97 If C{options} is not present, the options are reset to default 98 @param options: L{SolverOptions} 99 @note: if the adaption of subtolerance is choosen, the tolerance set by C{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 M{(Q^*Q)p=Q^*G} 107 108 @return: L{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 M{(Q^*Q)p=Q^*G} 116 117 If C{options} is not present, the options are reset to default 118 @param options: L{SolverOptions} 119 @note: if the adaption of subtolerance is choosen, the tolerance set by C{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. L{Data}) 129 @param g: flux sources/sinks 130 @type g: vector values on the domain (e.g. L{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. L{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. L{Data}) 135 @param permeability: permeability tensor. If scalar C{s} is given the tensor with 136 C{s} on the main diagonal is used. If vector C{v} is given the tensor with 137 C{v} on the main diagonal is used. 138 @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data}) 139 140 @note: the values of parameters which are not set by calling C{setValue} are not altered. 141 @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0) 142 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal 143 is along the M{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 C{rtol} used to terminate the solution process. The iteration is terminated if 181 182 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) } 183 184 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}}. 185 186 @param rtol: relative tolerance for the pressure 187 @type rtol: non-negative C{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: C{float} 198 """ 199 return self.__rtol 200 201 def setAbsoluteTolerance(self,atol=0.): 202 """ 203 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if 204 205 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) } 206 207 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}}. 208 209 @param atol: absolute tolerance for the pressure 210 @type atol: non-negative C{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: C{float} 221 """ 222 return self.__atol 223 def getSubProblemTolerance(self): 224 """ 225 Returns a suitable subtolerance 226 @type: C{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 C{location_of_fixed_flux} the value of C{u0} is kept unchanged. 248 @type u0: vector value on the domain (e.g. L{Data}). 249 @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. 250 @type p0: scalar value on the domain (e.g. L{Data}). 251 @param verbose: if set some information on iteration progress are printed 252 @type verbose: C{bool} 253 @return: flux and pressure 254 @rtype: C{tuple} of L{Data}. 255 256 @note: The problem is solved as a least squares form 257 258 M{(I+D^*D)u+Qp=D^*f+g} 259 M{Q^*u+Q^*Qp=Q^*g} 260 261 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. 262 We eliminate the flux form the problem by setting 263 264 M{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 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 269 270 which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step 271 PDEs with operator M{I+D^*D} and with M{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 278 num_corrections=0 279 converged=False 280 p=p0 281 norm_r=None 282 while not converged: 283 v=self.getFlux(p, fixed_flux=u0) 284 Qp=self.__Q(p) 285 norm_v=self.__L2(v) 286 norm_Qp=self.__L2(Qp) 287 if norm_v == 0.: 288 if norm_Qp == 0.: 289 return v,p 290 else: 291 fac=norm_Qp 292 else: 293 if norm_Qp == 0.: 294 fac=norm_v 295 else: 296 fac=2./(1./norm_v+1./norm_Qp) 297 ATOL=(atol+rtol*fac) 298 if self.verbose: 299 print "DarcyFlux: L2 norm of v = %e."%norm_v 300 print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp 301 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),) 302 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),) 303 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL 304 if norm_r == None or norm_r>ATOL: 305 if num_corrections>max_num_corrections: 306 raise ValueError,"maximum number of correction steps reached." 307 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) 308 num_corrections+=1 309 else: 310 converged=True 311 return v,p 312 def __L2(self,v): 313 return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2)) 314 315 def __Q(self,p): 316 return util.tensor_mult(self.__permeability,util.grad(p)) 317 318 def __Aprod(self,dp): 319 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator" 320 Qdp=self.__Q(dp) 321 self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data()) 322 du=self.__pde_v.getSolution() 323 # self.__pde_v.getOperator().saveMM("proj.mm") 324 return Qdp+du 325 def __inner_GMRES(self,r,s): 326 return util.integrate(util.inner(r,s)) 327 328 def __inner_PCG(self,p,r): 329 return util.integrate(util.inner(self.__Q(p), r)) 330 331 def __Msolve_PCG(self,r): 332 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner" 333 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data()) 334 # self.__pde_p.getOperator().saveMM("prec.mm") 335 return self.__pde_p.getSolution() 336 337 def getFlux(self,p=None, fixed_flux=Data()): 338 """ 339 returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux} 340 on locations where C{location_of_fixed_flux} is positive (see L{setValue}). 341 Note that C{g} and C{f} are used, see L{setValue}. 342 343 @param p: pressure. 344 @type p: scalar value on the domain (e.g. L{Data}). 345 @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}. 346 @type fixed_flux: vector values on the domain (e.g. L{Data}). 347 @param tol: relative tolerance to be used. 348 @type tol: positive C{float}. 349 @return: flux 350 @rtype: L{Data} 351 @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}} 352 for the permeability M{k_{ij}} 353 """ 354 self.setSubProblemTolerance() 355 g=self.__g 356 f=self.__f 357 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux) 358 if p == None: 359 self.__pde_v.setValue(Y=g) 360 else: 361 self.__pde_v.setValue(Y=g-self.__Q(p)) 362 return self.__pde_v.getSolution() 363 364 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 365 """ 366 solves 367 368 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j} 369 u_{i,i}=0 370 371 u=0 where fixed_u_mask>0 372 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 373 374 if surface_stress is not given 0 is assumed. 375 376 typical usage: 377 378 sp=StokesProblemCartesian(domain) 379 sp.setTolerance() 380 sp.initialize(...) 381 v,p=sp.solve(v0,p0) 382 """ 383 def __init__(self,domain,adaptSubTolerance=True, **kwargs): 384 """ 385 initialize the Stokes Problem 386 387 @param domain: domain of the problem. The approximation order needs to be two. 388 @type domain: L{Domain} 389 @param adaptSubTolerance: If True the tolerance for subproblem is set automatically. 390 @type adaptSubTolerance: C{bool} 391 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure. 392 """ 393 HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs) 394 self.domain=domain 395 self.vol=util.integrate(1.,Function(self.domain)) 396 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 397 self.__pde_u.setSymmetryOn() 398 399 self.__pde_prec=LinearPDE(domain) 400 self.__pde_prec.setReducedOrderOn() 401 self.__pde_prec.setSymmetryOn() 402 403 self.__pde_proj=LinearPDE(domain) 404 self.__pde_proj.setReducedOrderOn() 405 self.__pde_proj.setValue(D=1) 406 self.__pde_proj.setSymmetryOn() 407 408 def getSolverOptionsVelocity(self): 409 """ 410 returns the solver options used solve the equation for velocity. 411 412 @rtype: L{SolverOptions} 413 """ 414 return self.__pde_u.getSolverOptions() 415 def setSolverOptionsVelocity(self, options=None): 416 """ 417 set the solver options for solving the equation for velocity. 418 419 @param options: new solver options 420 @type options: L{SolverOptions} 421 """ 422 self.__pde_u.setSolverOptions(options) 423 def getSolverOptionsPressure(self): 424 """ 425 returns the solver options used solve the equation for pressure. 426 @rtype: L{SolverOptions} 427 """ 428 return self.__pde_prec.getSolverOptions() 429 def setSolverOptionsPressure(self, options=None): 430 """ 431 set the solver options for solving the equation for pressure. 432 @param options: new solver options 433 @type options: L{SolverOptions} 434 """ 435 self.__pde_prec.setSolverOptions(options) 436 437 def setSolverOptionsDiv(self, options=None): 438 """ 439 set the solver options for solving the equation to project the divergence of 440 the velocity onto the function space of presure. 441 442 @param options: new solver options 443 @type options: L{SolverOptions} 444 """ 445 self.__pde_prec.setSolverOptions(options) 446 def getSolverOptionsDiv(self): 447 """ 448 returns the solver options for solving the equation to project the divergence of 449 the velocity onto the function space of presure. 450 451 @rtype: L{SolverOptions} 452 """ 453 return self.__pde_prec.getSolverOptions() 454 def setSubProblemTolerance(self): 455 """ 456 Updates the tolerance for subproblems 457 """ 458 if self.adaptSubTolerance(): 459 sub_tol=self.getSubProblemTolerance() 460 self.getSolverOptionsDiv().setTolerance(sub_tol) 461 self.getSolverOptionsDiv().setAbsoluteTolerance(0.) 462 self.getSolverOptionsPressure().setTolerance(sub_tol) 463 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 464 self.getSolverOptionsVelocity().setTolerance(sub_tol) 465 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.) 466 467 468 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()): 469 """ 470 assigns values to the model parameters 471 472 @param f: external force 473 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar 474 @param fixed_u_mask: mask of locations with fixed velocity. 475 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar 476 @param eta: viscosity 477 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar 478 @param surface_stress: normal surface stress 479 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar 480 @param stress: initial stress 481 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar 482 @note: All values needs to be set. 483 """ 484 self.eta=eta 485 A =self.__pde_u.createCoefficient("A") 486 self.__pde_u.setValue(A=Data()) 487 for i in range(self.domain.getDim()): 488 for j in range(self.domain.getDim()): 489 A[i,j,j,i] += 1. 490 A[i,j,i,j] += 1. 491 self.__pde_prec.setValue(D=1/self.eta) 492 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask) 493 self.__f=f 494 self.__surface_stress=surface_stress 495 self.__stress=stress 496 497 def Bv(self,v): 498 """ 499 returns inner product of element p and div(v) 500 501 @param p: a pressure increment 502 @param v: a residual 503 @return: inner product of element p and div(v) 504 @rtype: C{float} 505 """ 506 self.__pde_proj.setValue(Y=-util.div(v)) 507 return self.__pde_proj.getSolution() 508 509 def inner_pBv(self,p,Bv): 510 """ 511 returns inner product of element p and Bv=-div(v) 512 513 @param p: a pressure increment 514 @param v: a residual 515 @return: inner product of element p and Bv=-div(v) 516 @rtype: C{float} 517 """ 518 return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain))) 519 520 def inner_p(self,p0,p1): 521 """ 522 Returns inner product of p0 and p1 523 524 @param p0: a pressure 525 @param p1: a pressure 526 @return: inner product of p0 and p1 527 @rtype: C{float} 528 """ 529 s0=util.interpolate(p0/self.eta,Function(self.domain)) 530 s1=util.interpolate(p1/self.eta,Function(self.domain)) 531 return util.integrate(s0*s1) 532 533 def norm_v(self,v): 534 """ 535 returns the norm of v 536 537 @param v: a velovity 538 @return: norm of v 539 @rtype: non-negative C{float} 540 """ 541 return util.sqrt(util.integrate(util.length(util.grad(v)))) 542 543 def getV(self, p, v0): 544 """ 545 return the value for v for a given p (overwrite) 546 547 @param p: a pressure 548 @param v0: a initial guess for the value v to return. 549 @return: v given as M{v= A^{-1} (f-B^*p)} 550 """ 551 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0) 552 if self.__stress.isEmpty(): 553 self.__pde_u.setValue(X=p*util.kronecker(self.domain)) 554 else: 555 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)) 556 out=self.__pde_u.getSolution() 557 return out 558 559 def norm_Bv(self,Bv): 560 """ 561 Returns Bv (overwrite). 562 563 @rtype: equal to the type of p 564 @note: boundary conditions on p should be zero! 565 """ 566 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2)) 567 568 def solve_AinvBt(self,p): 569 """ 570 Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()} 571 572 @param p: a pressure increment 573 @return: the solution of M{Av=B^*p} 574 @note: boundary conditions on v should be zero! 575 """ 576 self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain)) 577 out=self.__pde_u.getSolution() 578 return out 579 580 def solve_prec(self,Bv): 581 """ 582 applies preconditioner for for M{BA^{-1}B^*} to M{Bv} 583 with accuracy L{self.getSubProblemTolerance()} 584 585 @param v: velocity increment 586 @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*} 587 @note: boundary conditions on p are zero. 588 """ 589 self.__pde_prec.setValue(Y=Bv) 590 return self.__pde_prec.getSolution()

 ViewVC Help Powered by ViewVC 1.1.26