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

Revision 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 15980 byte(s)
```This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module.
Some documentation for both classes has been added.
The convection code is only linear at the moment.

```
 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 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): 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=boundary_flux) 141 return self.__pde_v.getSolution() 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: 202 raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, util.sqrt(norm_rhs), atol) 203 rhs=ArithmeticTuple(g,dv) 204 dp,r=PCG(rhs,self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=p0, verbose=self.verbose, initial_guess=False) 205 return u2+r[1],p2+dp 206 207 def __Aprod_PCG(self,p): 208 if self.show_details: print "DarcyFlux: Applying operator" 209 Qp=util.tensor_mult(self.__permeability,util.grad(p)) 210 self.__pde_v.setValue(Y=Qp,X=Data()) 211 w=self.__pde_v.getSolution(verbose=self.show_details) 212 return ArithmeticTuple(Qp,w) 213 214 def __inner_PCG(self,p,r): 215 a=util.tensor_mult(self.__permeability,util.grad(p)) 216 return util.integrate(util.inner(a,r[0]-r[1])) 217 218 def __Msolve_PCG(self,r): 219 if self.show_details: print "DarcyFlux: Applying preconditioner" 220 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1])) 221 return self.__pde_p.getSolution(verbose=self.show_details) 222 223 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 224 """ 225 solves 226 227 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j} 228 u_{i,i}=0 229 230 u=0 where fixed_u_mask>0 231 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 232 233 if surface_stress is not given 0 is assumed. 234 235 typical usage: 236 237 sp=StokesProblemCartesian(domain) 238 sp.setTolerance() 239 sp.initialize(...) 240 v,p=sp.solve(v0,p0) 241 """ 242 def __init__(self,domain,**kwargs): 243 """ 244 initialize the Stokes Problem 245 246 @param domain: domain of the problem. The approximation order needs to be two. 247 @type domain: L{Domain} 248 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure. 249 """ 250 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 251 self.domain=domain 252 self.vol=util.integrate(1.,Function(self.domain)) 253 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 254 self.__pde_u.setSymmetryOn() 255 # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT) 256 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU) 257 258 self.__pde_prec=LinearPDE(domain) 259 self.__pde_prec.setReducedOrderOn() 260 self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING) 261 self.__pde_prec.setSymmetryOn() 262 263 self.__pde_proj=LinearPDE(domain) 264 self.__pde_proj.setReducedOrderOn() 265 self.__pde_proj.setSymmetryOn() 266 self.__pde_proj.setValue(D=1.) 267 268 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()): 269 """ 270 assigns values to the model parameters 271 272 @param f: external force 273 @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar 274 @param fixed_u_mask: mask of locations with fixed velocity. 275 @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar 276 @param eta: viscosity 277 @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar 278 @param surface_stress: normal surface stress 279 @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar 280 @param stress: initial stress 281 @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar 282 @note: All values needs to be set. 283 284 """ 285 self.eta=eta 286 A =self.__pde_u.createCoefficient("A") 287 self.__pde_u.setValue(A=Data()) 288 for i in range(self.domain.getDim()): 289 for j in range(self.domain.getDim()): 290 A[i,j,j,i] += 1. 291 A[i,j,i,j] += 1. 292 self.__pde_prec.setValue(D=1/self.eta) 293 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) 294 self.__stress=stress 295 296 def B(self,v): 297 """ 298 returns div(v) 299 @rtype: equal to the type of p 300 301 @note: boundary conditions on p should be zero! 302 """ 303 if self.show_details: print "apply divergence:" 304 self.__pde_proj.setValue(Y=-util.div(v)) 305 self.__pde_proj.setTolerance(self.getSubProblemTolerance()) 306 return self.__pde_proj.getSolution(verbose=self.show_details) 307 308 def inner_pBv(self,p,Bv): 309 """ 310 returns inner product of element p and Bv (overwrite) 311 312 @type p: equal to the type of p 313 @type Bv: equal to the type of result of operator B 314 @rtype: C{float} 315 316 @rtype: equal to the type of p 317 """ 318 s0=util.interpolate(p,Function(self.domain)) 319 s1=util.interpolate(Bv,Function(self.domain)) 320 return util.integrate(s0*s1) 321 322 def inner_p(self,p0,p1): 323 """ 324 returns inner product of element p0 and p1 (overwrite) 325 326 @type p0: equal to the type of p 327 @type p1: equal to the type of p 328 @rtype: C{float} 329 330 @rtype: equal to the type of p 331 """ 332 s0=util.interpolate(p0/self.eta,Function(self.domain)) 333 s1=util.interpolate(p1/self.eta,Function(self.domain)) 334 return util.integrate(s0*s1) 335 336 def inner_v(self,v0,v1): 337 """ 338 returns inner product of two element v0 and v1 (overwrite) 339 340 @type v0: equal to the type of v 341 @type v1: equal to the type of v 342 @rtype: C{float} 343 344 @rtype: equal to the type of v 345 """ 346 gv0=util.grad(v0) 347 gv1=util.grad(v1) 348 return util.integrate(util.inner(gv0,gv1)) 349 350 def solve_A(self,u,p): 351 """ 352 solves Av=f-Au-B^*p (v=0 on fixed_u_mask) 353 """ 354 if self.show_details: print "solve for velocity:" 355 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 356 if self.__stress.isEmpty(): 357 self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain)) 358 else: 359 self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain)) 360 out=self.__pde_u.getSolution(verbose=self.show_details) 361 return out 362 363 def solve_prec(self,p): 364 if self.show_details: print "apply preconditioner:" 365 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 366 self.__pde_prec.setValue(Y=p) 367 q=self.__pde_prec.getSolution(verbose=self.show_details) 368 return q

 ViewVC Help Powered by ViewVC 1.1.26