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

Revision 2960 - (show annotations)
Tue Mar 2 07:54:11 2010 UTC (9 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 39649 byte(s)
```new DarcyFlux solver
```
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2010 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-2010 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 weight: the weighting factor for the incempressiblity equation *u_{i,i} = f* in the variational principle. 57 If not present an apprppriate weight is chosen. 58 :type weight: positive sacalar 59 :param useReduced: uses reduced oreder on flux and pressure 60 :type useReduced: ``bool`` 61 :param adaptSubTolerance: switches on automatic subtolerance selection 62 :type adaptSubTolerance: ``bool`` 63 """ 64 self.domain=domain 65 if weight == None: 66 self.__l=None 67 self.__update_weight=True 68 else: 69 self.__update_weight=False 70 self.__l=weight 71 self.__pde_v=LinearPDESystem(domain) 72 if useReduced: self.__pde_v.setReducedOrderOn() 73 self.__pde_v.setSymmetryOn() 74 self.__pde_p=LinearSinglePDE(domain) 75 self.__pde_p.setSymmetryOn() 76 if useReduced: self.__pde_p.setReducedOrderOn() 77 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 78 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 79 self.setTolerance() 80 self.setAbsoluteTolerance() 81 self.__adaptSubTolerance=adaptSubTolerance 82 self.verbose=False 83 def getSolverOptionsFlux(self): 84 """ 85 Returns the solver options used to solve the flux problems 86 87 *(I+D^* weight D)u=F* 88 89 :return: `SolverOptions` 90 """ 91 return self.__pde_v.getSolverOptions() 92 93 def setSolverOptionsFlux(self, options=None): 94 """ 95 Sets the solver options used to solve the flux problems 96 97 *(I+D^* weight D)u=F* 98 99 If ``options`` is not present, the options are reset to default 100 :param options: `SolverOptions` 101 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 102 """ 103 return self.__pde_v.setSolverOptions(options) 104 def getSolverOptionsPressure(self): 105 """ 106 Returns the solver options used to solve the pressure problems 107 108 *(Q^* S Q)p=Q^*G* 109 110 as a preconditioner where S=weight*k^2/dx^2 111 112 :return: `SolverOptions` 113 """ 114 return self.__pde_p.getSolverOptions() 115 def setSolverOptionsPressure(self, options=None): 116 """ 117 Sets the solver options used to solve the pressure problems 118 119 *(Q^* S Q)p=Q^*G* 120 121 was a preconditioner here S=weight*k^2/dx^2 122 123 If ``options`` is not present, the options are reset to default 124 125 :param options: `SolverOptions` 126 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 127 """ 128 return self.__pde_p.setSolverOptions(options) 129 130 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 131 """ 132 assigns values to model parameters 133 134 :param f: volumetic sources/sinks 135 :type f: scalar value on the domain (e.g. `Data`) 136 :param g: flux sources/sinks 137 :type g: vector values on the domain (e.g. `Data`) 138 :param location_of_fixed_pressure: mask for locations where pressure is fixed 139 :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`) 140 :param location_of_fixed_flux: mask for locations where flux is fixed. 141 :type location_of_fixed_flux: vector values on the domain (e.g. `Data`) 142 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with 143 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with 144 ``v`` on the main diagonal is used. 145 :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`) 146 147 :note: the values of parameters which are not set by calling ``setValue`` are not altered. 148 :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0) 149 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal 150 is along the *x_i* axis. 151 """ 152 if f !=None: 153 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X")) 154 if f.isEmpty(): 155 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 156 else: 157 if f.getRank()>0: raise ValueError,"illegal rank of f." 158 self.__f=f 159 if g !=None: 160 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y")) 161 if g.isEmpty(): 162 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 163 else: 164 if not g.getShape()==(self.domain.getDim(),): 165 raise ValueError,"illegal shape of g" 166 self.__g=g 167 168 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 169 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux) 170 171 if permeability!=None: 172 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A")) 173 if perm.getRank()==0: 174 perm_inv=1./perm*util.kronecker(self.domain.getDim()) 175 perm=perm*util.kronecker(self.domain.getDim()) 176 elif perm.getRank()==1: 177 perm2=perm 178 perm=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")) 179 perm_inv=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")) 180 for i in range(self.domain.getDim()): 181 p=perm2[i] 182 perm[i,i]=p 183 perm_inv[i,i]=1/p 184 elif perm.getRank()==2: 185 perm_inv=util.inverse(perm) 186 else: 187 raise ValueError,"illegal rank of permeability." 188 self.__permeability=perm 189 self.__permeability_inv=perm_inv 190 if self.__update_weight: self.__l = util.longestEdge(self.domain)**2/3.14**2*util.length(self.__permeability_inv) 191 self.__pde_v.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain))) 192 s=self.__pde_p.getFunctionSpaceForCoefficient("A").getSize() 193 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)*self.__l/s**2) 194 195 def getFlux(self,p=None, fixed_flux=Data()): 196 """ 197 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux`` 198 on locations where ``location_of_fixed_flux`` is positive (see `setValue`). 199 Note that ``g`` and ``f`` are used, see `setValue`. 200 201 :param p: pressure. 202 :type p: scalar value on the domain (e.g. `Data`). 203 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``. 204 :type fixed_flux: vector values on the domain (e.g. `Data`). 205 :return: flux 206 :rtype: `Data` 207 :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}* 208 for the permeability *k_{ij}* 209 """ 210 self.setSubProblemTolerance() 211 self.__pde_v.setValue(X=util.sqrt(self.__l)*self.__f*util.kronecker(self.domain), r=fixed_flux) 212 g2=util.tensor_mult(self.__permeability_inv,self.__g) 213 if p == None: 214 self.__pde_v.setValue(Y=g2) 215 else: 216 self.__pde_v.setValue(Y=g2-util.grad(p)) 217 return self.__pde_v.getSolution() 218 219 def __Aprod(self,dp): 220 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator" 221 Gdp=util.grad(dp) 222 self.__pde_v.setValue(Y=-Gdp,X=Data(), r=Data()) 223 du=self.__pde_v.getSolution() 224 # self.__pde_v.getOperator().saveMM("proj.mm") 225 return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du) 226 227 def __Msolve_PCG(self,r): 228 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner" 229 self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data()) 230 # self.__pde_p.getOperator().saveMM("prec.mm") 231 return self.__pde_p.getSolution() 232 233 def __inner_PCG(self,p,r): 234 return util.integrate(util.inner(util.grad(p), r[0]-r[1])) 235 236 def __L2(self,v): 237 return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2)) 238 239 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10): 240 """ 241 solves the problem. 242 243 The iteration is terminated if the residual norm is less then self.getTolerance(). 244 245 :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. 246 :type u0: vector value on the domain (e.g. `Data`). 247 :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. 248 :type p0: scalar value on the domain (e.g. `Data`). 249 :param verbose: if set some information on iteration progress are printed 250 :type verbose: ``bool`` 251 :return: flux and pressure 252 :rtype: ``tuple`` of `Data`. 253 254 :note: The problem is solved as a least squares form 255 256 *(K^[-1]+D^* weight D)u+G p=D^* sqrt(weight) f + K^[-1]g* 257 *G^*u+G^* K Gp=G^*g* 258 259 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*. 260 We eliminate the flux form the problem by setting 261 262 *u=(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g* -G p ) with u=u0 on location_of_fixed_flux 263 264 form the first equation. Inserted into the second equation we get 265 266 *G^*(K-(K^[-1]+D^* weight D)^{-1})Gp= G^*(g-(K^[-1]+D^* weight D)^{-1}(D^* sqrt(weight) f + K^[-1]g))* with p=p0 on location_of_fixed_pressure 267 268 which is solved using the PCG method (precondition is *G^*K^2*weight*l^2/dx^2 G*. In each iteration step 269 PDEs with operator *K^[-1]+D^* weight D* and with *G^*K^2*weight*l^2/dx^2 G* needs to be solved using a sub iteration scheme. 270 """ 271 self.verbose=verbose 272 rtol=self.getTolerance() 273 atol=self.getAbsoluteTolerance() 274 self.setSubProblemTolerance() 275 num_corrections=0 276 converged=False 277 p=p0 278 norm_r=None 279 v=self.getFlux(p, fixed_flux=u0) 280 while not converged: 281 Gp=util.grad(p) 282 KGp=util.tensor_mult(self.__permeability,Gp) 283 norm_v=util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5 284 norm_Gp=util.integrate(util.inner(Gp,KGp))**0.5 285 if self.verbose: 286 print "DarcyFlux: L2: g-v-K*grad(p) = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-KGp)**2)**(0.5),) 287 print "DarcyFlux: L2: f-div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),) 288 print "DarcyFlux: K^{-1}-norm of v = %e."%norm_v 289 print "DarcyFlux: K-norm of grad(p) = %e."%norm_Gp 290 if norm_v == 0.: 291 if norm_Gp == 0.: 292 return v,p 293 else: 294 fac=norm_Gp 295 else: 296 if norm_Gp == 0.: 297 fac=norm_v 298 else: 299 fac=2./(1./norm_v+1./norm_Gp) 300 #fac=max(norm_Gp, norm_v) 301 ATOL=atol+rtol*fac 302 if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e (fac= %e)."%(ATOL,fac) 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 # initial residual is r=G^*(self.__g-KGp - v) 307 p,r, norm_r=PCG(ArithmeticTuple(self.__g-KGp,v), 308 self.__Aprod, 309 p, 310 self.__Msolve_PCG, 311 self.__inner_PCG, 312 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose) 313 if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r 314 v=r[1] 315 num_corrections+=1 316 else: 317 if self.verbose: print "DarcyFlux: stopping criterium reached." 318 converged=True 319 return v,p 320 321 322 323 def setTolerance(self,rtol=1e-4): 324 """ 325 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if 326 327 *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}* 328 329 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`). 330 331 :param rtol: relative tolerance for the pressure 332 :type rtol: non-negative ``float`` 333 """ 334 if rtol<0: 335 raise ValueError,"Relative tolerance needs to be non-negative." 336 self.__rtol=rtol 337 def getTolerance(self): 338 """ 339 returns the relative tolerance 340 341 :return: current relative tolerance 342 :rtype: ``float`` 343 """ 344 return self.__rtol 345 346 def setAbsoluteTolerance(self,atol=0.): 347 """ 348 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if 349 350 *|g-v-K gard(p)|_PCG <= atol + rtol * (1/|K^{-1/2}u|_0 +1/|K^{1/2}grad(p)|_0)^(-1}* 351 352 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}*. 353 354 :param atol: absolute tolerance for the pressure 355 :type atol: non-negative ``float`` 356 """ 357 if atol<0: 358 raise ValueError,"Absolute tolerance needs to be non-negative." 359 self.__atol=atol 360 def getAbsoluteTolerance(self): 361 """ 362 returns the absolute tolerance 363 364 :return: current absolute tolerance 365 :rtype: ``float`` 366 """ 367 return self.__atol 368 def getSubProblemTolerance(self): 369 """ 370 Returns a suitable subtolerance 371 @type: ``float`` 372 """ 373 return max(util.EPSILON**(0.75),self.getTolerance()**2) 374 375 def setSubProblemTolerance(self): 376 """ 377 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected. 378 """ 379 if self.__adaptSubTolerance: 380 sub_tol=self.getSubProblemTolerance() 381 self.getSolverOptionsFlux().setTolerance(sub_tol) 382 self.getSolverOptionsFlux().setAbsoluteTolerance(0.) 383 self.getSolverOptionsPressure().setTolerance(sub_tol) 384 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 385 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol 386 387 388 class DarcyFlowOld(object): 389 """ 390 solves the problem 391 392 *u_i+k_{ij}*p_{,j} = g_i* 393 *u_{i,i} = f* 394 395 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability, 396 397 :note: The problem is solved in a least squares formulation. 398 """ 399 400 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True): 401 """ 402 initializes the Darcy flux problem 403 :param domain: domain of the problem 404 :type domain: `Domain` 405 :param useReduced: uses reduced oreder on flux and pressure 406 :type useReduced: ``bool`` 407 :param adaptSubTolerance: switches on automatic subtolerance selection 408 :type adaptSubTolerance: ``bool`` 409 """ 410 self.domain=domain 411 if weight == None: 412 s=self.domain.getSize() 413 self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2 414 # self.__l=(3.*util.longestEdge(self.domain))**2 415 #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2 416 else: 417 self.__l=weight 418 self.__pde_v=LinearPDESystem(domain) 419 if useReduced: self.__pde_v.setReducedOrderOn() 420 self.__pde_v.setSymmetryOn() 421 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain))) 422 self.__pde_p=LinearSinglePDE(domain) 423 self.__pde_p.setSymmetryOn() 424 if useReduced: self.__pde_p.setReducedOrderOn() 425 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 426 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 427 self.setTolerance() 428 self.setAbsoluteTolerance() 429 self.__adaptSubTolerance=adaptSubTolerance 430 self.verbose=False 431 def getSolverOptionsFlux(self): 432 """ 433 Returns the solver options used to solve the flux problems 434 435 *(I+D^*D)u=F* 436 437 :return: `SolverOptions` 438 """ 439 return self.__pde_v.getSolverOptions() 440 def setSolverOptionsFlux(self, options=None): 441 """ 442 Sets the solver options used to solve the flux problems 443 444 *(I+D^*D)u=F* 445 446 If ``options`` is not present, the options are reset to default 447 :param options: `SolverOptions` 448 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 449 """ 450 return self.__pde_v.setSolverOptions(options) 451 def getSolverOptionsPressure(self): 452 """ 453 Returns the solver options used to solve the pressure problems 454 455 *(Q^*Q)p=Q^*G* 456 457 :return: `SolverOptions` 458 """ 459 return self.__pde_p.getSolverOptions() 460 def setSolverOptionsPressure(self, options=None): 461 """ 462 Sets the solver options used to solve the pressure problems 463 464 *(Q^*Q)p=Q^*G* 465 466 If ``options`` is not present, the options are reset to default 467 :param options: `SolverOptions` 468 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 469 """ 470 return self.__pde_p.setSolverOptions(options) 471 472 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 473 """ 474 assigns values to model parameters 475 476 :param f: volumetic sources/sinks 477 :type f: scalar value on the domain (e.g. `Data`) 478 :param g: flux sources/sinks 479 :type g: vector values on the domain (e.g. `Data`) 480 :param location_of_fixed_pressure: mask for locations where pressure is fixed 481 :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`) 482 :param location_of_fixed_flux: mask for locations where flux is fixed. 483 :type location_of_fixed_flux: vector values on the domain (e.g. `Data`) 484 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with 485 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with 486 ``v`` on the main diagonal is used. 487 :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`) 488 489 :note: the values of parameters which are not set by calling ``setValue`` are not altered. 490 :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0) 491 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal 492 is along the *x_i* axis. 493 """ 494 if f !=None: 495 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X")) 496 if f.isEmpty(): 497 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 498 else: 499 if f.getRank()>0: raise ValueError,"illegal rank of f." 500 self.__f=f 501 if g !=None: 502 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y")) 503 if g.isEmpty(): 504 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 505 else: 506 if not g.getShape()==(self.domain.getDim(),): 507 raise ValueError,"illegal shape of g" 508 self.__g=g 509 510 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 511 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux) 512 513 if permeability!=None: 514 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A")) 515 if perm.getRank()==0: 516 perm=perm*util.kronecker(self.domain.getDim()) 517 elif perm.getRank()==1: 518 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm 519 for i in range(self.domain.getDim()): perm[i,i]=perm2[i] 520 elif perm.getRank()==2: 521 pass 522 else: 523 raise ValueError,"illegal rank of permeability." 524 self.__permeability=perm 525 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)) 526 527 def setTolerance(self,rtol=1e-4): 528 """ 529 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if 530 531 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )* 532 533 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}*. 534 535 :param rtol: relative tolerance for the pressure 536 :type rtol: non-negative ``float`` 537 """ 538 if rtol<0: 539 raise ValueError,"Relative tolerance needs to be non-negative." 540 self.__rtol=rtol 541 def getTolerance(self): 542 """ 543 returns the relative tolerance 544 545 :return: current relative tolerance 546 :rtype: ``float`` 547 """ 548 return self.__rtol 549 550 def setAbsoluteTolerance(self,atol=0.): 551 """ 552 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if 553 554 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )* 555 556 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}*. 557 558 :param atol: absolute tolerance for the pressure 559 :type atol: non-negative ``float`` 560 """ 561 if atol<0: 562 raise ValueError,"Absolute tolerance needs to be non-negative." 563 self.__atol=atol 564 def getAbsoluteTolerance(self): 565 """ 566 returns the absolute tolerance 567 568 :return: current absolute tolerance 569 :rtype: ``float`` 570 """ 571 return self.__atol 572 def getSubProblemTolerance(self): 573 """ 574 Returns a suitable subtolerance 575 @type: ``float`` 576 """ 577 return max(util.EPSILON**(0.75),self.getTolerance()**2) 578 def setSubProblemTolerance(self): 579 """ 580 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected. 581 """ 582 if self.__adaptSubTolerance: 583 sub_tol=self.getSubProblemTolerance() 584 self.getSolverOptionsFlux().setTolerance(sub_tol) 585 self.getSolverOptionsFlux().setAbsoluteTolerance(0.) 586 self.getSolverOptionsPressure().setTolerance(sub_tol) 587 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 588 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol 589 590 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10): 591 """ 592 solves the problem. 593 594 The iteration is terminated if the residual norm is less then self.getTolerance(). 595 596 :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. 597 :type u0: vector value on the domain (e.g. `Data`). 598 :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. 599 :type p0: scalar value on the domain (e.g. `Data`). 600 :param verbose: if set some information on iteration progress are printed 601 :type verbose: ``bool`` 602 :return: flux and pressure 603 :rtype: ``tuple`` of `Data`. 604 605 :note: The problem is solved as a least squares form 606 607 *(I+D^*D)u+Qp=D^*f+g* 608 *Q^*u+Q^*Qp=Q^*g* 609 610 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*. 611 We eliminate the flux form the problem by setting 612 613 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux 614 615 form the first equation. Inserted into the second equation we get 616 617 *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0 on location_of_fixed_pressure 618 619 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step 620 PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme. 621 """ 622 self.verbose=verbose 623 rtol=self.getTolerance() 624 atol=self.getAbsoluteTolerance() 625 self.setSubProblemTolerance() 626 num_corrections=0 627 converged=False 628 p=p0 629 norm_r=None 630 while not converged: 631 v=self.getFlux(p, fixed_flux=u0) 632 Qp=self.__Q(p) 633 norm_v=self.__L2(v) 634 norm_Qp=self.__L2(Qp) 635 if norm_v == 0.: 636 if norm_Qp == 0.: 637 return v,p 638 else: 639 fac=norm_Qp 640 else: 641 if norm_Qp == 0.: 642 fac=norm_v 643 else: 644 fac=2./(1./norm_v+1./norm_Qp) 645 ATOL=(atol+rtol*fac) 646 if self.verbose: 647 print "DarcyFlux: L2 norm of v = %e."%norm_v 648 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp 649 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),) 650 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),) 651 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL 652 if norm_r == None or norm_r>ATOL: 653 if num_corrections>max_num_corrections: 654 raise ValueError,"maximum number of correction steps reached." 655 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) 656 num_corrections+=1 657 else: 658 converged=True 659 return v,p 660 def __L2(self,v): 661 return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2)) 662 663 def __Q(self,p): 664 return util.tensor_mult(self.__permeability,util.grad(p)) 665 666 def __Aprod(self,dp): 667 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator" 668 Qdp=self.__Q(dp) 669 self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data()) 670 du=self.__pde_v.getSolution() 671 # self.__pde_v.getOperator().saveMM("proj.mm") 672 return Qdp+du 673 def __inner_GMRES(self,r,s): 674 return util.integrate(util.inner(r,s)) 675 676 def __inner_PCG(self,p,r): 677 return util.integrate(util.inner(self.__Q(p), r)) 678 679 def __Msolve_PCG(self,r): 680 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner" 681 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data()) 682 # self.__pde_p.getOperator().saveMM("prec.mm") 683 return self.__pde_p.getSolution() 684 685 def getFlux(self,p=None, fixed_flux=Data()): 686 """ 687 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux`` 688 on locations where ``location_of_fixed_flux`` is positive (see `setValue`). 689 Note that ``g`` and ``f`` are used, see `setValue`. 690 691 :param p: pressure. 692 :type p: scalar value on the domain (e.g. `Data`). 693 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``. 694 :type fixed_flux: vector values on the domain (e.g. `Data`). 695 :return: flux 696 :rtype: `Data` 697 :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}* 698 for the permeability *k_{ij}* 699 """ 700 self.setSubProblemTolerance() 701 g=self.__g 702 f=self.__f 703 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux) 704 if p == None: 705 self.__pde_v.setValue(Y=g) 706 else: 707 self.__pde_v.setValue(Y=g-self.__Q(p)) 708 return self.__pde_v.getSolution() 709 710 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 711 """ 712 solves 713 714 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j} 715 u_{i,i}=0 716 717 u=0 where fixed_u_mask>0 718 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 719 720 if surface_stress is not given 0 is assumed. 721 722 typical usage: 723 724 sp=StokesProblemCartesian(domain) 725 sp.setTolerance() 726 sp.initialize(...) 727 v,p=sp.solve(v0,p0) 728 """ 729 def __init__(self,domain,**kwargs): 730 """ 731 initialize the Stokes Problem 732 733 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be 734 LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation 735 with macro elements for the pressure. 736 737 :param domain: domain of the problem. 738 :type domain: `Domain` 739 """ 740 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 741 self.domain=domain 742 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 743 self.__pde_u.setSymmetryOn() 744 745 self.__pde_prec=LinearPDE(domain) 746 self.__pde_prec.setReducedOrderOn() 747 self.__pde_prec.setSymmetryOn() 748 749 self.__pde_proj=LinearPDE(domain) 750 self.__pde_proj.setReducedOrderOn() 751 self.__pde_proj.setValue(D=1) 752 self.__pde_proj.setSymmetryOn() 753 754 def getSolverOptionsVelocity(self): 755 """ 756 returns the solver options used solve the equation for velocity. 757 758 :rtype: `SolverOptions` 759 """ 760 return self.__pde_u.getSolverOptions() 761 def setSolverOptionsVelocity(self, options=None): 762 """ 763 set the solver options for solving the equation for velocity. 764 765 :param options: new solver options 766 :type options: `SolverOptions` 767 """ 768 self.__pde_u.setSolverOptions(options) 769 def getSolverOptionsPressure(self): 770 """ 771 returns the solver options used solve the equation for pressure. 772 :rtype: `SolverOptions` 773 """ 774 return self.__pde_prec.getSolverOptions() 775 def setSolverOptionsPressure(self, options=None): 776 """ 777 set the solver options for solving the equation for pressure. 778 :param options: new solver options 779 :type options: `SolverOptions` 780 """ 781 self.__pde_prec.setSolverOptions(options) 782 783 def setSolverOptionsDiv(self, options=None): 784 """ 785 set the solver options for solving the equation to project the divergence of 786 the velocity onto the function space of presure. 787 788 :param options: new solver options 789 :type options: `SolverOptions` 790 """ 791 self.__pde_proj.setSolverOptions(options) 792 def getSolverOptionsDiv(self): 793 """ 794 returns the solver options for solving the equation to project the divergence of 795 the velocity onto the function space of presure. 796 797 :rtype: `SolverOptions` 798 """ 799 return self.__pde_proj.getSolverOptions() 800 801 def updateStokesEquation(self, v, p): 802 """ 803 updates the Stokes equation to consider dependencies from ``v`` and ``p`` 804 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values. 805 """ 806 pass 807 def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None): 808 """ 809 assigns new values to the model parameters. 810 811 :param f: external force 812 :type f: `Vector` object in `FunctionSpace` `Function` or similar 813 :param fixed_u_mask: mask of locations with fixed velocity. 814 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar 815 :param eta: viscosity 816 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar 817 :param surface_stress: normal surface stress 818 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar 819 :param stress: initial stress 820 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar 821 """ 822 if eta !=None: 823 k=util.kronecker(self.domain.getDim()) 824 kk=util.outer(k,k) 825 self.eta=util.interpolate(eta, Function(self.domain)) 826 self.__pde_prec.setValue(D=1/self.eta) 827 self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3))) 828 if restoration_factor!=None: 829 n=self.domain.getNormal() 830 self.__pde_u.setValue(d=restoration_factor*util.outer(n,n)) 831 if fixed_u_mask!=None: 832 self.__pde_u.setValue(q=fixed_u_mask) 833 if f!=None: self.__f=f 834 if surface_stress!=None: self.__surface_stress=surface_stress 835 if stress!=None: self.__stress=stress 836 837 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0): 838 """ 839 assigns values to the model parameters 840 841 :param f: external force 842 :type f: `Vector` object in `FunctionSpace` `Function` or similar 843 :param fixed_u_mask: mask of locations with fixed velocity. 844 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar 845 :param eta: viscosity 846 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar 847 :param surface_stress: normal surface stress 848 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar 849 :param stress: initial stress 850 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar 851 """ 852 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor) 853 854 def Bv(self,v,tol): 855 """ 856 returns inner product of element p and div(v) 857 858 :param v: a residual 859 :return: inner product of element p and div(v) 860 :rtype: ``float`` 861 """ 862 self.__pde_proj.setValue(Y=-util.div(v)) 863 self.getSolverOptionsDiv().setTolerance(tol) 864 self.getSolverOptionsDiv().setAbsoluteTolerance(0.) 865 out=self.__pde_proj.getSolution() 866 return out 867 868 def inner_pBv(self,p,Bv): 869 """ 870 returns inner product of element p and Bv=-div(v) 871 872 :param p: a pressure increment 873 :param Bv: a residual 874 :return: inner product of element p and Bv=-div(v) 875 :rtype: ``float`` 876 """ 877 return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain))) 878 879 def inner_p(self,p0,p1): 880 """ 881 Returns inner product of p0 and p1 882 883 :param p0: a pressure 884 :param p1: a pressure 885 :return: inner product of p0 and p1 886 :rtype: ``float`` 887 """ 888 s0=util.interpolate(p0,Function(self.domain)) 889 s1=util.interpolate(p1,Function(self.domain)) 890 return util.integrate(s0*s1) 891 892 def norm_v(self,v): 893 """ 894 returns the norm of v 895 896 :param v: a velovity 897 :return: norm of v 898 :rtype: non-negative ``float`` 899 """ 900 return util.sqrt(util.integrate(util.length(util.grad(v))**2)) 901 902 903 def getDV(self, p, v, tol): 904 """ 905 return the value for v for a given p (overwrite) 906 907 :param p: a pressure 908 :param v: a initial guess for the value v to return. 909 :return: dv given as *Adv=(f-Av-B^*p)* 910 """ 911 self.updateStokesEquation(v,p) 912 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress) 913 self.getSolverOptionsVelocity().setTolerance(tol) 914 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.) 915 if self.__stress.isEmpty(): 916 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 917 else: 918 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 919 out=self.__pde_u.getSolution() 920 return out 921 922 def norm_Bv(self,Bv): 923 """ 924 Returns Bv (overwrite). 925 926 :rtype: equal to the type of p 927 :note: boundary conditions on p should be zero! 928 """ 929 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2)) 930 931 def solve_AinvBt(self,p, tol): 932 """ 933 Solves *Av=B^*p* with accuracy `tol` 934 935 :param p: a pressure increment 936 :return: the solution of *Av=B^*p* 937 :note: boundary conditions on v should be zero! 938 """ 939 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain)) 940 out=self.__pde_u.getSolution() 941 return out 942 943 def solve_prec(self,Bv, tol): 944 """ 945 applies preconditioner for for *BA^{-1}B^** to *Bv* 946 with accuracy `self.getSubProblemTolerance()` 947 948 :param Bv: velocity increment 949 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )* 950 :note: boundary conditions on p are zero. 951 """ 952 self.__pde_prec.setValue(Y=Bv) 953 self.getSolverOptionsPressure().setTolerance(tol) 954 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 955 out=self.__pde_prec.getSolution() 956 return out

 ViewVC Help Powered by ViewVC 1.1.26