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

Revision 2156 - (show annotations)
Mon Dec 15 05:09:02 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 17615 byte(s)
```some modifications to the iterative solver to make them easier to use.
There are also improved versions of the Darcy flux solver and the incompressible solver.

```
 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 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): 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.__pde_v=LinearPDESystem(domain) 59 self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain))) 60 self.__pde_v.setSymmetryOn() 61 self.__pde_p=LinearSinglePDE(domain) 62 self.__pde_p.setSymmetryOn() 63 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 64 self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 65 66 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 67 """ 68 assigns values to model parameters 69 70 @param f: volumetic sources/sinks 71 @type f: scalar value on the domain (e.g. L{Data}) 72 @param g: flux sources/sinks 73 @type g: vector values on the domain (e.g. L{Data}) 74 @param location_of_fixed_pressure: mask for locations where pressure is fixed 75 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data}) 76 @param location_of_fixed_flux: mask for locations where flux is fixed. 77 @type location_of_fixed_flux: vector values on the domain (e.g. L{Data}) 78 @param permeability: permeability tensor. If scalar C{s} is given the tensor with 79 C{s} on the main diagonal is used. If vector C{v} is given the tensor with 80 C{v} on the main diagonal is used. 81 @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data}) 82 83 @note: the values of parameters which are not set by calling C{setValue} are not altered. 84 @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0) 85 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal 86 is along the M{x_i} axis. 87 """ 88 if f !=None: 89 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X")) 90 if f.isEmpty(): 91 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 92 else: 93 if f.getRank()>0: raise ValueError,"illegal rank of f." 94 self.f=f 95 if g !=None: 96 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y")) 97 if g.isEmpty(): 98 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 99 else: 100 if not g.getShape()==(self.domain.getDim(),): 101 raise ValueError,"illegal shape of g" 102 self.__g=g 103 104 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 105 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux) 106 107 if permeability!=None: 108 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A")) 109 if perm.getRank()==0: 110 perm=perm*util.kronecker(self.domain.getDim()) 111 elif perm.getRank()==1: 112 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm 113 for i in range(self.domain.getDim()): perm[i,i]=perm2[i] 114 elif perm.getRank()==2: 115 pass 116 else: 117 raise ValueError,"illegal rank of permeability." 118 self.__permeability=perm 119 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)) 120 121 122 def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False): 123 """ 124 returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux} 125 on locations where C{location_of_fixed_flux} is positive (see L{setValue}). 126 Note that C{g} and C{f} are used, L{setValue}. 127 128 @param p: pressure. 129 @type p: scalar value on the domain (e.g. L{Data}). 130 @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}. 131 @type fixed_flux: vector values on the domain (e.g. L{Data}). 132 @param tol: relative tolerance to be used. 133 @type tol: positive float. 134 @return: flux 135 @rtype: L{Data} 136 @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}} 137 for the permeability M{k_{ij}} 138 """ 139 self.__pde_v.setTolerance(tol) 140 self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux) 141 return self.__pde_v.getSolution(verbose=show_details) 142 143 def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8): 144 """ 145 solves the problem. 146 147 The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where 148 C{|q|} denotes the norm of the right hand side (see escript user's guide for details). 149 150 @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. 151 @type u0: vector value on the domain (e.g. L{Data}). 152 @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. 153 @type p0: scalar value on the domain (e.g. L{Data}). 154 @param atol: absolute tolerance for the pressure 155 @type atol: non-negative C{float} 156 @param rtol: relative tolerance for the pressure 157 @type rtol: non-negative C{float} 158 @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol0 265 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 266 267 if surface_stress is not given 0 is assumed. 268 269 typical usage: 270 271 sp=StokesProblemCartesian(domain) 272 sp.setTolerance() 273 sp.initialize(...) 274 v,p=sp.solve(v0,p0) 275 """ 276 def __init__(self,domain,**kwargs): 277 """ 278 initialize the Stokes Problem 279 280 @param domain: domain of the problem. The approximation order needs to be two. 281 @type domain: L{Domain} 282 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure. 283 """ 284 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 285 self.domain=domain 286 self.vol=util.integrate(1.,Function(self.domain)) 287 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 288 self.__pde_u.setSymmetryOn() 289 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT) 290 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU) 291 292 self.__pde_prec=LinearPDE(domain) 293 self.__pde_prec.setReducedOrderOn() 294 # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING) 295 self.__pde_prec.setSymmetryOn() 296 297 self.__pde_proj=LinearPDE(domain) 298 self.__pde_proj.setReducedOrderOn() 299 self.__pde_proj.setSymmetryOn() 300 self.__pde_proj.setValue(D=1.) 301 302 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()): 303 """ 304 assigns values to the model parameters 305 306 @param f: external force 307 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar 308 @param fixed_u_mask: mask of locations with fixed velocity. 309 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar 310 @param eta: viscosity 311 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar 312 @param surface_stress: normal surface stress 313 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar 314 @param stress: initial stress 315 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar 316 @note: All values needs to be set. 317 318 """ 319 self.eta=eta 320 A =self.__pde_u.createCoefficient("A") 321 self.__pde_u.setValue(A=Data()) 322 for i in range(self.domain.getDim()): 323 for j in range(self.domain.getDim()): 324 A[i,j,j,i] += 1. 325 A[i,j,i,j] += 1. 326 self.__pde_prec.setValue(D=1/self.eta) 327 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) 328 self.__stress=stress 329 330 def B(self,v): 331 """ 332 returns div(v) 333 @rtype: equal to the type of p 334 335 @note: boundary conditions on p should be zero! 336 """ 337 if self.show_details: print "apply divergence:" 338 self.__pde_proj.setValue(Y=-util.div(v)) 339 self.__pde_proj.setTolerance(self.getSubProblemTolerance()) 340 return self.__pde_proj.getSolution(verbose=self.show_details) 341 342 def inner_pBv(self,p,Bv): 343 """ 344 returns inner product of element p and Bv (overwrite) 345 346 @type p: equal to the type of p 347 @type Bv: equal to the type of result of operator B 348 @rtype: C{float} 349 350 @rtype: equal to the type of p 351 """ 352 s0=util.interpolate(p,Function(self.domain)) 353 s1=util.interpolate(Bv,Function(self.domain)) 354 return util.integrate(s0*s1) 355 356 def inner_p(self,p0,p1): 357 """ 358 returns inner product of element p0 and p1 (overwrite) 359 360 @type p0: equal to the type of p 361 @type p1: equal to the type of p 362 @rtype: C{float} 363 364 @rtype: equal to the type of p 365 """ 366 s0=util.interpolate(p0/self.eta,Function(self.domain)) 367 s1=util.interpolate(p1/self.eta,Function(self.domain)) 368 return util.integrate(s0*s1) 369 370 def inner_v(self,v0,v1): 371 """ 372 returns inner product of two element v0 and v1 (overwrite) 373 374 @type v0: equal to the type of v 375 @type v1: equal to the type of v 376 @rtype: C{float} 377 378 @rtype: equal to the type of v 379 """ 380 gv0=util.grad(v0) 381 gv1=util.grad(v1) 382 return util.integrate(util.inner(gv0,gv1)) 383 384 def solve_A(self,u,p): 385 """ 386 solves Av=f-Au-B^*p (v=0 on fixed_u_mask) 387 """ 388 if self.show_details: print "solve for velocity:" 389 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 390 if self.__stress.isEmpty(): 391 self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain)) 392 else: 393 self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain)) 394 out=self.__pde_u.getSolution(verbose=self.show_details) 395 return out 396 397 def solve_prec(self,p): 398 if self.show_details: print "apply preconditioner:" 399 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 400 self.__pde_prec.setValue(Y=p) 401 q=self.__pde_prec.getSolution(verbose=self.show_details) 402 return q

 ViewVC Help Powered by ViewVC 1.1.26