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

Revision 2719 - (show annotations)
Wed Oct 14 06:38:03 2009 UTC (11 years ago) by gross
File MIME type: text/x-python
File size: 22832 byte(s)
```a new Stokes solver added
```
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2009 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-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 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 *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,**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 :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure. 387 """ 388 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 389 self.domain=domain 390 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 391 self.__pde_u.setSymmetryOn() 392 393 self.__pde_prec=LinearPDE(domain) 394 self.__pde_prec.setReducedOrderOn() 395 self.__pde_prec.setSymmetryOn() 396 397 self.__pde_proj=LinearPDE(domain) 398 self.__pde_proj.setReducedOrderOn() 399 self.__pde_proj.setValue(D=1) 400 self.__pde_proj.setSymmetryOn() 401 402 def getSolverOptionsVelocity(self): 403 """ 404 returns the solver options used solve the equation for velocity. 405 406 :rtype: `SolverOptions` 407 """ 408 return self.__pde_u.getSolverOptions() 409 def setSolverOptionsVelocity(self, options=None): 410 """ 411 set the solver options for solving the equation for velocity. 412 413 :param options: new solver options 414 :type options: `SolverOptions` 415 """ 416 self.__pde_u.setSolverOptions(options) 417 def getSolverOptionsPressure(self): 418 """ 419 returns the solver options used solve the equation for pressure. 420 :rtype: `SolverOptions` 421 """ 422 return self.__pde_prec.getSolverOptions() 423 def setSolverOptionsPressure(self, options=None): 424 """ 425 set the solver options for solving the equation for pressure. 426 :param options: new solver options 427 :type options: `SolverOptions` 428 """ 429 self.__pde_prec.setSolverOptions(options) 430 431 def setSolverOptionsDiv(self, options=None): 432 """ 433 set the solver options for solving the equation to project the divergence of 434 the velocity onto the function space of presure. 435 436 :param options: new solver options 437 :type options: `SolverOptions` 438 """ 439 self.__pde_proj.setSolverOptions(options) 440 def getSolverOptionsDiv(self): 441 """ 442 returns the solver options for solving the equation to project the divergence of 443 the velocity onto the function space of presure. 444 445 :rtype: `SolverOptions` 446 """ 447 return self.__pde_proj.getSolverOptions() 448 449 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0): 450 """ 451 assigns values to the model parameters 452 453 :param f: external force 454 :type f: `Vector` object in `FunctionSpace` `Function` or similar 455 :param fixed_u_mask: mask of locations with fixed velocity. 456 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar 457 :param eta: viscosity 458 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar 459 :param surface_stress: normal surface stress 460 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar 461 :param stress: initial stress 462 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar 463 :note: All values needs to be set. 464 """ 465 self.eta=eta 466 A =self.__pde_u.createCoefficient("A") 467 self.__pde_u.setValue(A=Data()) 468 for i in range(self.domain.getDim()): 469 for j in range(self.domain.getDim()): 470 A[i,j,j,i] += 1. 471 A[i,j,i,j] += 1. 472 n=self.domain.getNormal() 473 self.__pde_prec.setValue(D=1/self.eta) 474 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n)) 475 self.__f=f 476 self.__surface_stress=surface_stress 477 self.__stress=stress 478 479 def Bv(self,v,tol): 480 """ 481 returns inner product of element p and div(v) 482 483 :param v: a residual 484 :return: inner product of element p and div(v) 485 :rtype: ``float`` 486 """ 487 self.__pde_proj.setValue(Y=-util.div(v)) # -??? 488 self.getSolverOptionsDiv().setTolerance(tol) 489 self.getSolverOptionsDiv().setAbsoluteTolerance(0.) 490 out=self.__pde_proj.getSolution() 491 return out 492 493 def inner_pBv(self,p,Bv): 494 """ 495 returns inner product of element p and Bv=-div(v) 496 497 :param p: a pressure increment 498 :param Bv: a residual 499 :return: inner product of element p and Bv=-div(v) 500 :rtype: ``float`` 501 """ 502 return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain))) 503 504 def inner_p(self,p0,p1): 505 """ 506 Returns inner product of p0 and p1 507 508 :param p0: a pressure 509 :param p1: a pressure 510 :return: inner product of p0 and p1 511 :rtype: ``float`` 512 """ 513 s0=util.interpolate(p0,Function(self.domain)) 514 s1=util.interpolate(p1,Function(self.domain)) 515 return util.integrate(s0*s1) 516 517 def norm_v(self,v): 518 """ 519 returns the norm of v 520 521 :param v: a velovity 522 :return: norm of v 523 :rtype: non-negative ``float`` 524 """ 525 return util.sqrt(util.integrate(util.length(util.grad(v))**2)) 526 527 def getDV(self, p, v, tol): 528 """ 529 return the value for v for a given p (overwrite) 530 531 :param p: a pressure 532 :param v: a initial guess for the value v to return. 533 :return: dv given as *Adv=(f-Av-B^*p)* 534 """ 535 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress) 536 self.getSolverOptionsVelocity().setTolerance(tol) 537 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.) 538 if self.__stress.isEmpty(): 539 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 540 else: 541 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 542 out=self.__pde_u.getSolution() 543 return out 544 545 def norm_Bv(self,Bv): 546 """ 547 Returns Bv (overwrite). 548 549 :rtype: equal to the type of p 550 :note: boundary conditions on p should be zero! 551 """ 552 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2)) 553 554 def solve_AinvBt(self,p, tol): 555 """ 556 Solves *Av=B^*p* with accuracy `tol` 557 558 :param p: a pressure increment 559 :return: the solution of *Av=B^*p* 560 :note: boundary conditions on v should be zero! 561 """ 562 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain)) 563 out=self.__pde_u.getSolution() 564 return out 565 566 def solve_prec(self,Bv, tol): 567 """ 568 applies preconditioner for for *BA^{-1}B^** to *Bv* 569 with accuracy `self.getSubProblemTolerance()` 570 571 :param Bv: velocity increment 572 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )* 573 :note: boundary conditions on p are zero. 574 """ 575 self.__pde_prec.setValue(Y=Bv) 576 self.getSolverOptionsPressure().setTolerance(tol) 577 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 578 out=self.__pde_prec.getSolution() 579 return out

 ViewVC Help Powered by ViewVC 1.1.26