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

Revision 2267 - (show annotations)
Thu Feb 12 06:18:45 2009 UTC (10 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 22467 byte(s)
weighting factor introduced.
 1 ######################################################## 2 # 3 # Copyright (c) 2003-2008 by University of Queensland 4 # Earth Systems Science Computational Center (ESSCC) 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 ######################################################## 12 13 __copyright__="""Copyright (c) 2003-2008 by University of Queensland 14 Earth Systems Science Computational Center (ESSCC) 15 http://www.uq.edu.au/esscc 16 Primary Business: Queensland, Australia""" 17 __license__="""Licensed under the Open Software License version 3.0 18 19 __url__= 20 21 """ 22 Some models for flow 23 24 @var __author__: name of author 25 @var __copyright__: copyrights 26 @var __license__: licence agreement 27 @var __url__: url entry point on documentation 28 @var __version__: version 29 @var __date__: date of the version 30 """ 31 32 __author__="Lutz Gross, l.gross@uq.edu.au" 33 34 from escript import * 35 import util 36 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE 37 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES 38 39 class DarcyFlow(object): 40 """ 41 solves the problem 42 43 M{u_i+k_{ij}*p_{,j} = g_i} 44 M{u_{i,i} = f} 45 46 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability, 47 48 @note: The problem is solved in a least squares formulation. 49 """ 50 51 def __init__(self, domain,useReduced=False): 52 """ 53 initializes the Darcy flux problem 54 @param domain: domain of the problem 55 @type domain: L{Domain} 56 """ 57 self.domain=domain 58 self.__l=util.integrate(Scalar(1.,Function(self.domain)))*(1./self.domain.getDim()) 59 self.__pde_v=LinearPDESystem(domain) 60 if useReduced: self.__pde_v.setReducedOrderOn() 61 self.__pde_v.setSymmetryOn() 62 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain))) 63 self.__pde_p=LinearSinglePDE(domain) 64 self.__pde_p.setSymmetryOn() 65 if useReduced: self.__pde_p.setReducedOrderOn() 66 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 67 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 68 self.setTolerance() 69 self.setAbsoluteTolerance() 70 self.setSubProblemTolerance() 71 72 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 73 """ 74 assigns values to model parameters 75 76 @param f: volumetic sources/sinks 77 @type f: scalar value on the domain (e.g. L{Data}) 78 @param g: flux sources/sinks 79 @type g: vector values on the domain (e.g. L{Data}) 80 @param location_of_fixed_pressure: mask for locations where pressure is fixed 81 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data}) 82 @param location_of_fixed_flux: mask for locations where flux is fixed. 83 @type location_of_fixed_flux: vector values on the domain (e.g. L{Data}) 84 @param permeability: permeability tensor. If scalar C{s} is given the tensor with 85 C{s} on the main diagonal is used. If vector C{v} is given the tensor with 86 C{v} on the main diagonal is used. 87 @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data}) 88 89 @note: the values of parameters which are not set by calling C{setValue} are not altered. 90 @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0) 91 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal 92 is along the M{x_i} axis. 93 """ 94 if f !=None: 95 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X")) 96 if f.isEmpty(): 97 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 98 else: 99 if f.getRank()>0: raise ValueError,"illegal rank of f." 100 self.__f=f 101 if g !=None: 102 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y")) 103 if g.isEmpty(): 104 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 105 else: 106 if not g.getShape()==(self.domain.getDim(),): 107 raise ValueError,"illegal shape of g" 108 self.__g=g 109 110 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 111 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux) 112 113 if permeability!=None: 114 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A")) 115 if perm.getRank()==0: 116 perm=perm*util.kronecker(self.domain.getDim()) 117 elif perm.getRank()==1: 118 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm 119 for i in range(self.domain.getDim()): perm[i,i]=perm2[i] 120 elif perm.getRank()==2: 121 pass 122 else: 123 raise ValueError,"illegal rank of permeability." 124 self.__permeability=perm 125 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)) 126 127 def setTolerance(self,rtol=1e-4): 128 """ 129 sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if 130 131 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) } 132 133 where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. 134 135 @param rtol: relative tolerance for the pressure 136 @type rtol: non-negative C{float} 137 """ 138 if rtol<0: 139 raise ValueError,"Relative tolerance needs to be non-negative." 140 self.__rtol=rtol 141 def getTolerance(self): 142 """ 143 returns the relative tolerance 144 145 @return: current relative tolerance 146 @rtype: C{float} 147 """ 148 return self.__rtol 149 150 def setAbsoluteTolerance(self,atol=0.): 151 """ 152 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if 153 154 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) } 155 156 where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. 157 158 @param atol: absolute tolerance for the pressure 159 @type atol: non-negative C{float} 160 """ 161 if atol<0: 162 raise ValueError,"Absolute tolerance needs to be non-negative." 163 self.__atol=atol 164 def getAbsoluteTolerance(self): 165 """ 166 returns the absolute tolerance 167 168 @return: current absolute tolerance 169 @rtype: C{float} 170 """ 171 return self.__atol 172 173 def setSubProblemTolerance(self,rtol=None): 174 """ 175 Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present 176 C{self.getTolerance()**2} is used. 177 178 @param rtol: relative tolerence 179 @type rtol: positive C{float} 180 """ 181 if rtol == None: 182 if self.getTolerance()<=0.: 183 raise ValueError,"A positive relative tolerance must be set." 184 self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2) 185 else: 186 if rtol<=0: 187 raise ValueError,"sub-problem tolerance must be positive." 188 self.__sub_tol=max(util.EPSILON**(0.75),rtol) 189 190 def getSubProblemTolerance(self): 191 """ 192 Returns the subproblem reduction factor. 193 194 @return: subproblem reduction factor 195 @rtype: C{float} 196 """ 197 return self.__sub_tol 198 199 def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10): 200 """ 201 solves the problem. 202 203 The iteration is terminated if the residual norm is less then self.getTolerance(). 204 205 @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged. 206 @type u0: vector value on the domain (e.g. L{Data}). 207 @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged. 208 @type p0: scalar value on the domain (e.g. L{Data}). 209 @param verbose: if set some information on iteration progress are printed 210 @type verbose: C{bool} 211 @param show_details: if set information on the subiteration process are printed. 212 @type show_details: C{bool} 213 @return: flux and pressure 214 @rtype: C{tuple} of L{Data}. 215 216 @note: The problem is solved as a least squares form 217 218 M{(I+D^*D)u+Qp=D^*f+g} 219 M{Q^*u+Q^*Qp=Q^*g} 220 221 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. 222 We eliminate the flux form the problem by setting 223 224 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux 225 226 form the first equation. Inserted into the second equation we get 227 228 M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0 on location_of_fixed_pressure 229 230 which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step 231 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme. 232 """ 233 self.verbose=verbose 234 self.show_details= show_details and self.verbose 235 rtol=self.getTolerance() 236 atol=self.getAbsoluteTolerance() 237 if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance() 238 239 num_corrections=0 240 converged=False 241 p=p0 242 norm_r=None 243 while not converged: 244 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details) 245 Qp=self.__Q(p) 246 norm_v=self.__L2(v) 247 norm_Qp=self.__L2(Qp) 248 if norm_v == 0.: 249 if norm_Qp == 0.: 250 return v,p 251 else: 252 fac=norm_Qp 253 else: 254 if norm_Qp == 0.: 255 fac=norm_v 256 else: 257 fac=2./(1./norm_v+1./norm_Qp) 258 ATOL=(atol+rtol*fac) 259 if self.verbose: 260 print "DarcyFlux: L2 norm of v = %e."%norm_v 261 print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp 262 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL 263 if norm_r == None or norm_r>ATOL: 264 if num_corrections>max_num_corrections: 265 raise ValueError,"maximum number of correction steps reached." 266 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.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose) 267 num_corrections+=1 268 else: 269 converged=True 270 return v,p 271 # 272 # 273 # r_hat=g-util.interpolate(v,Function(self.domain))-Qp 274 # #=========================================================================== 275 # norm_r_hat=self.__L2(r_hat) 276 # norm_v=self.__L2(v) 277 # norm_g=self.__L2(g) 278 # norm_gv=self.__L2(g-v) 279 # norm_Qp=self.__L2(Qp) 280 # norm_gQp=self.__L2(g-Qp) 281 # fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv)) 282 # fac=min(norm_v,norm_Qp,norm_gv) 283 # norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat)) 284 # print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat 285 # if r!=None: 286 # print "diff = ",self.__L2(r-r_hat)/norm_r_hat 287 # sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance() 288 # self.setSubProblemTolerance(sub_tol) 289 # print "subtol_new=",self.getSubProblemTolerance() 290 # print "norm_v = ",norm_v 291 # print "norm_gv = ",norm_gv 292 # print "norm_Qp = ",norm_Qp 293 # print "norm_gQp = ",norm_gQp 294 # print "norm_g = ",norm_g 295 # print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp) 296 # print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv) 297 # if fac == 0: 298 # if self.verbose: print "DarcyFlux: trivial case!" 299 # return v,p 300 # #=============================================================================== 301 # # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v)) 302 # # norm_Qp=self.__L2(Qp) 303 # norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat)) 304 # # print "**** norm_v, norm_Qp :",norm_v,norm_Qp 305 # 306 # ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp)) 307 # if self.verbose: 308 # print "DarcyFlux: residual = %e"%norm_r_hat 309 # print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL 310 # if norm_r_hat <= ATOL: 311 # print "DarcyFlux: iteration finalized." 312 # converged=True 313 # else: 314 # # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG) 315 # # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose) 316 # p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose) 317 # print "norm_r =",norm_r 318 # return v,p 319 def __L2(self,v): 320 return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2)) 321 322 def __Q(self,p): 323 return util.tensor_mult(self.__permeability,util.grad(p)) 324 325 def __Aprod(self,dp): 326 self.__pde_v.setTolerance(self.getSubProblemTolerance()) 327 if self.show_details: print "DarcyFlux: Applying operator" 328 Qdp=self.__Q(dp) 329 self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data()) 330 du=self.__pde_v.getSolution(verbose=self.show_details) 331 return Qdp+du 332 def __inner_GMRES(self,r,s): 333 return util.integrate(util.inner(r,s)) 334 335 def __inner_PCG(self,p,r): 336 return util.integrate(util.inner(self.__Q(p), r)) 337 338 def __Msolve_PCG(self,r): 339 self.__pde_p.setTolerance(self.getSubProblemTolerance()) 340 if self.show_details: print "DarcyFlux: Applying preconditioner" 341 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data()) 342 return self.__pde_p.getSolution(verbose=self.show_details) 343 344 def getFlux(self,p=None, fixed_flux=Data(), show_details=False): 345 """ 346 returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux} 347 on locations where C{location_of_fixed_flux} is positive (see L{setValue}). 348 Note that C{g} and C{f} are used, see L{setValue}. 349 350 @param p: pressure. 351 @type p: scalar value on the domain (e.g. L{Data}). 352 @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}. 353 @type fixed_flux: vector values on the domain (e.g. L{Data}). 354 @param tol: relative tolerance to be used. 355 @type tol: positive C{float}. 356 @return: flux 357 @rtype: L{Data} 358 @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} 359 for the permeability M{k_{ij}} 360 """ 361 self.__pde_v.setTolerance(self.getSubProblemTolerance()) 362 g=self.__g 363 f=self.__f 364 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux) 365 if p == None: 366 self.__pde_v.setValue(Y=g) 367 else: 368 self.__pde_v.setValue(Y=g-self.__Q(p)) 369 return self.__pde_v.getSolution(verbose=show_details) 370 371 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 372 """ 373 solves 374 375 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j} 376 u_{i,i}=0 377 378 u=0 where fixed_u_mask>0 379 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 380 381 if surface_stress is not given 0 is assumed. 382 383 typical usage: 384 385 sp=StokesProblemCartesian(domain) 386 sp.setTolerance() 387 sp.initialize(...) 388 v,p=sp.solve(v0,p0) 389 """ 390 def __init__(self,domain,**kwargs): 391 """ 392 initialize the Stokes Problem 393 394 @param domain: domain of the problem. The approximation order needs to be two. 395 @type domain: L{Domain} 396 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure. 397 """ 398 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 399 self.domain=domain 400 self.vol=util.integrate(1.,Function(self.domain)) 401 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 402 self.__pde_u.setSymmetryOn() 403 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT) 404 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU) 405 406 self.__pde_prec=LinearPDE(domain) 407 self.__pde_prec.setReducedOrderOn() 408 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING) 409 self.__pde_prec.setSymmetryOn() 410 411 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()): 412 """ 413 assigns values to the model parameters 414 415 @param f: external force 416 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar 417 @param fixed_u_mask: mask of locations with fixed velocity. 418 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar 419 @param eta: viscosity 420 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar 421 @param surface_stress: normal surface stress 422 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar 423 @param stress: initial stress 424 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar 425 @note: All values needs to be set. 426 427 """ 428 self.eta=eta 429 A =self.__pde_u.createCoefficient("A") 430 self.__pde_u.setValue(A=Data()) 431 for i in range(self.domain.getDim()): 432 for j in range(self.domain.getDim()): 433 A[i,j,j,i] += 1. 434 A[i,j,i,j] += 1. 435 self.__pde_prec.setValue(D=1/self.eta) 436 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask) 437 self.__f=f 438 self.__surface_stress=surface_stress 439 self.__stress=stress 440 441 def inner_pBv(self,p,v): 442 """ 443 returns inner product of element p and div(v) 444 445 @param p: a pressure increment 446 @param v: a residual 447 @return: inner product of element p and div(v) 448 @rtype: C{float} 449 """ 450 return util.integrate(-p*util.div(v)) 451 452 def inner_p(self,p0,p1): 453 """ 454 Returns inner product of p0 and p1 455 456 @param p0: a pressure 457 @param p1: a pressure 458 @return: inner product of p0 and p1 459 @rtype: C{float} 460 """ 461 s0=util.interpolate(p0/self.eta,Function(self.domain)) 462 s1=util.interpolate(p1/self.eta,Function(self.domain)) 463 return util.integrate(s0*s1) 464 465 def norm_v(self,v): 466 """ 467 returns the norm of v 468 469 @param v: a velovity 470 @return: norm of v 471 @rtype: non-negative C{float} 472 """ 473 return util.sqrt(util.integrate(util.length(util.grad(v)))) 474 475 def getV(self, p, v0): 476 """ 477 return the value for v for a given p (overwrite) 478 479 @param p: a pressure 480 @param v0: a initial guess for the value v to return. 481 @return: v given as M{v= A^{-1} (f-B^*p)} 482 """ 483 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 484 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0) 485 if self.__stress.isEmpty(): 486 self.__pde_u.setValue(X=p*util.kronecker(self.domain)) 487 else: 488 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)) 489 out=self.__pde_u.getSolution(verbose=self.show_details) 490 return out 491 492 493 raise NotImplementedError,"no v calculation implemented." 494 495 496 def norm_Bv(self,v): 497 """ 498 Returns Bv (overwrite). 499 500 @rtype: equal to the type of p 501 @note: boundary conditions on p should be zero! 502 """ 503 return util.sqrt(util.integrate(util.div(v)**2)) 504 505 def solve_AinvBt(self,p): 506 """ 507 Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()} 508 509 @param p: a pressure increment 510 @return: the solution of M{Av=B^*p} 511 @note: boundary conditions on v should be zero! 512 """ 513 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 514 self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain)) 515 out=self.__pde_u.getSolution(verbose=self.show_details) 516 return out 517 518 def solve_precB(self,v): 519 """ 520 applies preconditioner for for M{BA^{-1}B^*} to M{Bv} 521 with accuracy L{self.getSubProblemTolerance()} (overwrite). 522 523 @param v: velocity increment 524 @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*} 525 @note: boundary conditions on p are zero. 526 """ 527 self.__pde_prec.setValue(Y=-util.div(v)) 528 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 529 return self.__pde_prec.getSolution(verbose=self.show_details)

 ViewVC Help Powered by ViewVC 1.1.26