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

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

 ViewVC Help Powered by ViewVC 1.1.26