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

Revision 3499 - (show annotations)
Fri Apr 8 00:14:35 2011 UTC (8 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 54972 byte(s)

 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 import escript 36 import util 37 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions 38 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES 39 40 41 class DarcyFlowVeryNew(object): 42 """ 43 solves the problem 44 45 *u_i+k_{ij}*p_{,j} = g_i* 46 *u_{i,i} = f* 47 48 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability, 49 50 :note: The problem is solved in a stabelized formulation. 51 """ 52 def __init__(self, domain, useReduced=False, *args, **kargs): 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 :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE. 64 :type useVPIteration: ``bool`` 65 """ 66 self.domain=domain 67 useVPIteration=False 68 self.useVPIteration=useVPIteration or True 69 #self.useVPIteration=False 70 self.useReduced=useReduced 71 self.verbose=False 72 73 if self.useVPIteration: 74 # this does not work yet 75 self.__pde_k=LinearPDESystem(domain) 76 self.__pde_k.setSymmetryOn() 77 if self.useReduced: self.__pde_k.setReducedOrderOn() 78 79 self.__pde_p=LinearSinglePDE(domain) 80 self.__pde_p.setSymmetryOn() 81 if self.useReduced: self.__pde_p.setReducedOrderOn() 82 else: 83 self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1) 84 self.__pde_k.setSymmetryOff() 85 86 if self.useReduced: self.__pde_k.setReducedOrderOn() 87 C=self.__pde_k.createCoefficient("C") 88 B=self.__pde_k.createCoefficient("B") 89 for i in range(self.domain.getDim()): 90 B[i,i,self.domain.getDim()]=-1 91 C[self.domain.getDim(),i,i]=1 92 C[i,self.domain.getDim(),i]=-0.5 93 B[self.domain.getDim(),i,i]=0.5 94 self.__pde_k.setValue(C=C, B=B) 95 self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X")) 96 self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y")) 97 self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q")) 98 self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q")) 99 self.scale=1. 100 def __L2(self,v): 101 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2)) 102 def __inner_GMRES(self,r,s): 103 return util.integrate(util.inner(r,s)) 104 105 def __Aprod_GMRES(self,p): 106 self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) ) 107 du=self.__pde_k.getSolution() 108 self.__pde_p.setValue(Y=util.div(du), X=0.5*(du+util.tensor_mult(self.__permeability,util.grad(p)))) 109 return self.__pde_p.getSolution() 110 111 def getSolverOptionsFlux(self): 112 """ 113 Returns the solver options used to solve the flux problems 114 115 *K^{-1} u=F* 116 117 :return: `SolverOptions` 118 """ 119 return self.__pde_k.getSolverOptions() 120 121 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 122 """ 123 assigns values to model parameters 124 125 :param f: volumetic sources/sinks 126 :type f: scalar value on the domain (e.g. `escript.Data`) 127 :param g: flux sources/sinks 128 :type g: vector values on the domain (e.g. `escript.Data`) 129 :param location_of_fixed_pressure: mask for locations where pressure is fixed 130 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`) 131 :param location_of_fixed_flux: mask for locations where flux is fixed. 132 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`) 133 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used. 134 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`) 135 136 :note: the values of parameters which are not set by calling ``setValue`` are not altered. 137 :note: at any point on the boundary of the domain the pressure 138 (``location_of_fixed_pressure`` >0) or the normal component of the 139 flux (``location_of_fixed_flux[i]>0``) if direction of the normal 140 is along the *x_i* axis. 141 142 """ 143 if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure 144 if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux 145 146 if self.useVPIteration: 147 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure) 148 if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux) 149 else: 150 if location_of_fixed_pressure!=None or location_of_fixed_flux!=None: 151 q=self.__pde_k.createCoefficient("q") 152 q[self.domain.getDim()]=self.location_of_fixed_pressure 153 q[:self.domain.getDim()]=self.location_of_fixed_flux 154 self.__pde_k.setValue(q=q) 155 156 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain) 157 if permeability!=None: 158 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A")) 159 V=util.vol(self.domain) 160 if perm.getRank()==0: 161 perm_inv=(1./perm) 162 self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.) 163 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim())) 164 perm=perm*(self.scale*util.kronecker(self.domain.getDim())) 165 elif perm.getRank()==2: 166 perm_inv=util.inverse(perm) 167 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim()) 168 perm_inv*=(1./self.scale) 169 perm=perm*self.scale 170 else: 171 raise ValueError,"illegal rank of permeability." 172 173 self.__permeability=perm 174 self.__permeability_inv=perm_inv 175 if self.useVPIteration: 176 self.__pde_k.setValue(D=0.5*self.__permeability_inv) 177 self.__pde_p.setValue(A=0.5*self.__permeability) 178 else: 179 D=self.__pde_k.createCoefficient("D") 180 A=self.__pde_k.createCoefficient("A") 181 D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv 182 A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability 183 self.__pde_k.setValue(A=A, D=D) 184 if g != None: 185 g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y")) 186 if g.isEmpty(): 187 g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y")) 188 else: 189 if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g" 190 self.__g=g 191 if f !=None: 192 f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X")) 193 if f.isEmpty(): 194 195 f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X")) 196 else: 197 if f.getRank()>0: raise ValueError,"illegal rank of f." 198 self.__f=f 199 200 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10): 201 """ 202 solves the problem. 203 204 The iteration is terminated if the residual norm is less then self.getTolerance(). 205 206 :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. 207 :type u0: vector value on the domain (e.g. `escript.Data`). 208 :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. 209 :type p0: scalar value on the domain (e.g. `escript.Data`). 210 :param verbose: if set some information on iteration progress are printed 211 :type verbose: ``bool`` 212 :return: flux and pressure 213 :rtype: ``tuple`` of `escript.Data`. 214 215 """ 216 u0_b=u0*self.location_of_fixed_flux 217 p0_b=p0*self.location_of_fixed_pressure/self.scale 218 f=self.__f-util.div(u0_b) 219 g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b)) 220 self.verbose=verbose 221 if self.useVPIteration: 222 # get u: 223 self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data()) 224 du=self.__pde_k.getSolution() 225 self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du)) 226 p=GMRES(self.__pde_p.getSolution(), 227 self.__Aprod_GMRES, 228 p0_b*0, 229 self.__inner_GMRES, 230 atol=0, 231 rtol=1.e-4, 232 iter_max=100, 233 iter_restart=20, verbose=self.verbose,P_R=None) 234 self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) , 235 X=p*util.kronecker(self.__pde_k.getDomain())) 236 u=self.__pde_k.getSolution() 237 else: 238 X=self.__pde_k.createCoefficient("X") 239 Y=self.__pde_k.createCoefficient("Y") 240 Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g) 241 Y[self.domain.getDim()]=f 242 X[self.domain.getDim(),:]=g*0.5 243 self.__pde_k.setValue(X=X, Y=Y) 244 self.__pde_k.getSolverOptions().setVerbosity(self.verbose) 245 #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG) 246 self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT) 247 U=self.__pde_k.getSolution() 248 u=U[:self.domain.getDim()] 249 p=U[self.domain.getDim()] 250 # self.__pde_k.getOperator().saveMM("k.mm") 251 u=u0_b+u 252 p=(p0_b+p)*self.scale 253 if self.verbose: 254 KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale) 255 def_p=self.__g-(u+KGp) 256 def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X")) 257 print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u)) 258 print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u))) 259 return u,p 260 def setTolerance(self,rtol=1e-4): 261 """ 262 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if 263 264 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0* 265 266 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`). 267 268 :param rtol: relative tolerance for the pressure 269 :type rtol: non-negative ``float`` 270 """ 271 if rtol<0: 272 raise ValueError,"Relative tolerance needs to be non-negative." 273 self.__rtol=rtol 274 def getTolerance(self): 275 """ 276 returns the relative tolerance 277 :return: current relative tolerance 278 :rtype: ``float`` 279 """ 280 return self.__rtol 281 class DarcyFlow(object): 282 """ 283 solves the problem 284 285 *u_i+k_{ij}*p_{,j} = g_i* 286 *u_{i,i} = f* 287 288 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability, 289 290 :note: The problem is solved in a least squares formulation. 291 """ 292 293 def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1): 294 """ 295 initializes the Darcy flux problem 296 :param domain: domain of the problem 297 :type domain: `Domain` 298 :param useReduced: uses reduced oreder on flux and pressure 299 :type useReduced: ``bool`` 300 :param adaptSubTolerance: switches on automatic subtolerance selection 301 :type adaptSubTolerance: ``bool`` 302 :param solveForFlux: if True the solver solves for the flux (do not use!) 303 :type solveForFlux: ``bool`` 304 :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE. 305 :type useVPIteration: ``bool`` 306 """ 307 self.domain=domain 308 self.useVPIteration=useVPIteration 309 self.useReduced=useReduced 310 self.weighting_scale=weighting_scale 311 if self.useVPIteration: 312 self.solveForFlux=solveForFlux 313 self.__adaptSubTolerance=adaptSubTolerance 314 self.verbose=False 315 316 self.__pde_k=LinearPDESystem(domain) 317 self.__pde_k.setSymmetryOn() 318 if self.useReduced: self.__pde_k.setReducedOrderOn() 319 320 self.__pde_p=LinearSinglePDE(domain) 321 self.__pde_p.setSymmetryOn() 322 if self.useReduced: self.__pde_p.setReducedOrderOn() 323 self.setTolerance() 324 self.setAbsoluteTolerance() 325 else: 326 self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1) 327 self.__pde_k.setSymmetryOn() 328 if self.useReduced: self.__pde_k.setReducedOrderOn() 329 C=self.__pde_k.createCoefficient("C") 330 B=self.__pde_k.createCoefficient("B") 331 for i in range(self.domain.getDim()): 332 C[i,self.domain.getDim(),i]=1 333 B[self.domain.getDim(),i,i]=1 334 self.__pde_k.setValue(C=C, B=B) 335 self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X")) 336 self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y")) 337 338 def getSolverOptionsFlux(self): 339 """ 340 Returns the solver options used to solve the flux problems 341 342 *K^{-1} u=F* 343 344 :return: `SolverOptions` 345 """ 346 return self.__pde_k.getSolverOptions() 347 348 def setSolverOptionsFlux(self, options=None): 349 """ 350 Sets the solver options used to solve the flux problems 351 352 *K^{-1}u=F* 353 354 If ``options`` is not present, the options are reset to default 355 356 :param options: `SolverOptions` 357 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 358 """ 359 return self.__pde_v.setSolverOptions(options) 360 361 def getSolverOptionsPressure(self): 362 """ 363 Returns the solver options used to solve the pressure problems 364 365 *(Q^* K Q)p=-Q^*G* 366 367 :return: `SolverOptions` 368 """ 369 return self.__pde_p.getSolverOptions() 370 371 def setSolverOptionsPressure(self, options=None): 372 """ 373 Sets the solver options used to solve the pressure problems 374 375 *(Q^* K Q)p=-Q^*G* 376 377 If ``options`` is not present, the options are reset to default 378 379 :param options: `SolverOptions` 380 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 381 """ 382 return self.__pde_p.setSolverOptions(options) 383 384 385 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 386 """ 387 assigns values to model parameters 388 389 :param f: volumetic sources/sinks 390 :type f: scalar value on the domain (e.g. `escript.Data`) 391 :param g: flux sources/sinks 392 :type g: vector values on the domain (e.g. `escript.Data`) 393 :param location_of_fixed_pressure: mask for locations where pressure is fixed 394 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`) 395 :param location_of_fixed_flux: mask for locations where flux is fixed. 396 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`) 397 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used. 398 :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`) 399 400 :note: the values of parameters which are not set by calling ``setValue`` are not altered. 401 :note: at any point on the boundary of the domain the pressure 402 (``location_of_fixed_pressure`` >0) or the normal component of the 403 flux (``location_of_fixed_flux[i]>0``) if direction of the normal 404 is along the *x_i* axis. 405 406 """ 407 if self.useVPIteration: 408 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 409 if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux) 410 else: 411 q=self.__pde_k.getCoefficient("q") 412 if q.isEmpty(): q=self.__pde_k.createCoefficient("q") 413 if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure 414 if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux 415 self.__pde_k.setValue(q=q) 416 417 # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain) 418 if permeability!=None: 419 perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A")) 420 V=util.vol(self.domain) 421 if perm.getRank()==0: 422 perm_inv=(1./perm) 423 if self.useVPIteration: 424 self.scale=1. 425 else: 426 self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.) 427 428 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim())) 429 perm=perm*(self.scale*util.kronecker(self.domain.getDim())) 430 elif perm.getRank()==2: 431 perm_inv=util.inverse(perm) 432 if self.useVPIteration: 433 self.scale=1. 434 else: 435 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim()) 436 perm_inv*=(1./self.scale) 437 perm=perm*self.scale 438 else: 439 raise ValueError,"illegal rank of permeability." 440 441 self.__permeability=perm 442 self.__permeability_inv=perm_inv 443 444 self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale 445 if self.useVPIteration: 446 if self.solveForFlux: 447 self.__pde_k.setValue(D=self.__permeability_inv) 448 else: 449 self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain))) 450 self.__pde_p.setValue(A=self.__permeability) 451 else: 452 D=self.__pde_k.createCoefficient("D") 453 A=self.__pde_k.createCoefficient("A") 454 D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv 455 for i in range(self.domain.getDim()): 456 for j in range(self.domain.getDim()): 457 A[i,i,j,j]=self.__l2 458 A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability 459 self.__pde_k.setValue(A=A, D=D) 460 if g !=None: 461 g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y")) 462 if g.isEmpty(): 463 g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y")) 464 else: 465 if not g.getShape()==(self.domain.getDim(),): 466 raise ValueError,"illegal shape of g" 467 self.__g=g 468 elif permeability!=None: 469 X 470 if f !=None: 471 f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X")) 472 if f.isEmpty(): 473 f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X")) 474 else: 475 if f.getRank()>0: raise ValueError,"illegal rank of f." 476 self.__f=f 477 478 479 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10): 480 """ 481 solves the problem. 482 483 The iteration is terminated if the residual norm is less then self.getTolerance(). 484 485 :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. 486 :type u0: vector value on the domain (e.g. `escript.Data`). 487 :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. 488 :type p0: scalar value on the domain (e.g. `escript.Data`). 489 :param verbose: if set some information on iteration progress are printed 490 :type verbose: ``bool`` 491 :return: flux and pressure 492 :rtype: ``tuple`` of `escript.Data`. 493 494 :note: The problem is solved as a least squares form 495 *(K^[-1]+D^* l2 D)u+G p=D^* l2 * f + K^[-1]g* 496 *G^*u+*G^* K Gp=G^*g* 497 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*. 498 """ 499 self.verbose=verbose 500 if self.useVPIteration: 501 return self.__solveVP(u0,p0,max_iter,max_num_corrections) 502 else: 503 X=self.__pde_k.createCoefficient("X") 504 Y=self.__pde_k.createCoefficient("Y") 505 Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g) 506 rtmp=self.__f * self.__l2 * self.scale 507 for i in range(self.domain.getDim()): X[i,i]=rtmp 508 X[self.domain.getDim(),:]=self.__g*self.scale 509 r=self.__pde_k.createCoefficient("r") 510 r[:self.domain.getDim()]=u0*self.scale 511 r[self.domain.getDim()]=p0 512 self.__pde_k.setValue(X=X, Y=Y, r=r) 513 self.__pde_k.getSolverOptions().setVerbosity(self.verbose) 514 #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG) 515 self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT) 516 U=self.__pde_k.getSolution() 517 # self.__pde_k.getOperator().saveMM("k.mm") 518 u=U[:self.domain.getDim()]*(1./self.scale) 519 p=U[self.domain.getDim()] 520 if self.verbose: 521 KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale) 522 def_p=self.__g-(u+KGp) 523 def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X")) 524 print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u)) 525 print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u))) 526 return u,p 527 528 def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10): 529 """ 530 solves the problem. 531 532 The iteration is terminated if the residual norm is less than self.getTolerance(). 533 534 :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. 535 :type u0: vector value on the domain (e.g. `escript.Data`). 536 :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. 537 :type p0: scalar value on the domain (e.g. `escript.Data`). 538 :return: flux and pressure 539 :rtype: ``tuple`` of `escript.Data`. 540 541 :note: The problem is solved as a least squares form 542 *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g* 543 *G^*u+*G^* K Gp=G^*g* 544 where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*. 545 """ 546 rtol=self.getTolerance() 547 atol=self.getAbsoluteTolerance() 548 self.setSubProblemTolerance() 549 num_corrections=0 550 converged=False 551 norm_r=None 552 553 # Eliminate the hydrostatic pressure: 554 if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component." 555 self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0)) 556 p0=self.__pde_p.getSolution() 557 g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0)) 558 norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5 559 560 p=p0*0 561 if self.solveForFlux: 562 v=u0.copy() 563 else: 564 v=self.__getFlux(p, u0, f=self.__f, g=g2) 565 566 while not converged and norm_g2 > 0: 567 Gp=util.grad(p) 568 KGp=util.tensor_mult(self.__permeability,Gp) 569 if self.verbose: 570 def_p=g2-(v+KGp) 571 def_v=self.__f-util.div(v) 572 print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v)) 573 print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v))) 574 print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5 575 print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2 576 print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5 577 ATOL=atol+rtol*norm_g2 578 if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,) 579 if norm_r == None or norm_r>ATOL: 580 if num_corrections>max_num_corrections: 581 raise ValueError,"maximum number of correction steps reached." 582 583 if self.solveForFlux: 584 # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du) 585 v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p), 586 self.__Aprod_v, 587 v, 588 self.__Msolve_PCG_v, 589 self.__inner_PCG_v, 590 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose) 591 p=r[2] 592 else: 593 # initial residual is r=G^*(g2-KGp - v) 594 p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v), 595 self.__Aprod_p, 596 p, 597 self.__Msolve_PCG_p, 598 self.__inner_PCG_p, 599 atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose) 600 v=r[1] 601 if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r 602 num_corrections+=1 603 else: 604 if self.verbose: print "DarcyFlux: stopping criterium reached." 605 converged=True 606 return v,p+p0 607 608 def __applWeight(self, v, f=None): 609 # solves L p = f-Dv with p = 0 610 if self.verbose: print "DarcyFlux: Applying weighting operator" 611 if f == None: 612 return -util.div(v)*self.__l2 613 else: 614 return (f-util.div(v))*self.__l2 615 def __getPressure(self, v, p0, g=None): 616 # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0 617 if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update" 618 if g == None: 619 self.__pde_p.setValue(X=-v, r=p0) 620 else: 621 self.__pde_p.setValue(X=g-v, r=p0) 622 p=self.__pde_p.getSolution() 623 return p 624 625 def __Aprod_v(self,dv): 626 # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp) with (G*KG)dp = - G^*dv 627 dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv) 628 a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp 629 b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv) 630 return ArithmeticTuple(a,b,-dp) 631 632 def __Msolve_PCG_v(self,r): 633 # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv 634 if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner" 635 self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=escript.Data()) 636 return self.__pde_k.getSolution() 637 638 def __inner_PCG_v(self,v,r): 639 return util.integrate(util.inner(v,r[0])+util.div(v)*r[1]) 640 641 def __Aprod_p(self,dp): 642 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator" 643 Gdp=util.grad(dp) 644 self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data()) 645 du=self.__pde_k.getSolution() 646 return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du) 647 648 def __getFlux(self,p, v0, f=None, g=None): 649 # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp 650 if f!=None: 651 self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain)) 652 self.__pde_k.setValue(r=v0) 653 g2=util.tensor_mult(self.__permeability_inv,g) 654 if p == None: 655 self.__pde_k.setValue(Y=g2) 656 else: 657 self.__pde_k.setValue(Y=g2-util.grad(p)) 658 return self.__pde_k.getSolution() 659 660 #v=self.__getFlux(p, u0, f=self.__f, g=g2) 661 def __Msolve_PCG_p(self,r): 662 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner" 663 self.__pde_p.setValue(X=r[0]-r[1], Y=escript.Data(), r=escript.Data(), y=escript.Data()) 664 return self.__pde_p.getSolution() 665 666 def __inner_PCG_p(self,p,r): 667 return util.integrate(util.inner(util.grad(p), r[0]-r[1])) 668 669 def __L2(self,v): 670 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2)) 671 672 def __L2_r(self,v): 673 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2)) 674 675 def setTolerance(self,rtol=1e-4): 676 """ 677 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if 678 679 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0* 680 681 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`). 682 683 :param rtol: relative tolerance for the pressure 684 :type rtol: non-negative ``float`` 685 """ 686 if rtol<0: 687 raise ValueError,"Relative tolerance needs to be non-negative." 688 self.__rtol=rtol 689 def getTolerance(self): 690 """ 691 returns the relative tolerance 692 :return: current relative tolerance 693 :rtype: ``float`` 694 """ 695 return self.__rtol 696 697 def setAbsoluteTolerance(self,atol=0.): 698 """ 699 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if 700 701 *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0* 702 703 704 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}*. 705 706 :param atol: absolute tolerance for the pressure 707 :type atol: non-negative ``float`` 708 """ 709 if atol<0: 710 raise ValueError,"Absolute tolerance needs to be non-negative." 711 self.__atol=atol 712 def getAbsoluteTolerance(self): 713 """ 714 returns the absolute tolerance 715 :return: current absolute tolerance 716 :rtype: ``float`` 717 """ 718 return self.__atol 719 def getSubProblemTolerance(self): 720 """ 721 Returns a suitable subtolerance 722 :type: ``float`` 723 """ 724 return max(util.EPSILON**(0.5),self.getTolerance()**2) 725 726 def setSubProblemTolerance(self): 727 """ 728 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected. 729 """ 730 if self.__adaptSubTolerance: 731 sub_tol=self.getSubProblemTolerance() 732 self.getSolverOptionsFlux().setTolerance(sub_tol) 733 self.getSolverOptionsFlux().setAbsoluteTolerance(0.) 734 self.getSolverOptionsPressure().setTolerance(sub_tol) 735 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 736 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol 737 738 739 class DarcyFlowOld(object): 740 """ 741 solves the problem 742 743 *u_i+k_{ij}*p_{,j} = g_i* 744 *u_{i,i} = f* 745 746 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability, 747 748 :note: The problem is solved in a least squares formulation. 749 """ 750 751 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True): 752 """ 753 initializes the Darcy flux problem 754 :param domain: domain of the problem 755 :type domain: `Domain` 756 :param useReduced: uses reduced oreder on flux and pressure 757 :type useReduced: ``bool`` 758 :param adaptSubTolerance: switches on automatic subtolerance selection 759 :type adaptSubTolerance: ``bool`` 760 """ 761 self.domain=domain 762 if weight == None: 763 s=self.domain.getSize() 764 self.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2 765 # self.__l2=(3.*util.longestEdge(self.domain))**2 766 #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2 767 else: 768 self.__l2=weight 769 self.__pde_v=LinearPDESystem(domain) 770 if useReduced: self.__pde_v.setReducedOrderOn() 771 self.__pde_v.setSymmetryOn() 772 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l2*util.outer(util.kronecker(domain),util.kronecker(domain))) 773 self.__pde_p=LinearSinglePDE(domain) 774 self.__pde_p.setSymmetryOn() 775 if useReduced: self.__pde_p.setReducedOrderOn() 776 self.__f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 777 self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 778 self.setTolerance() 779 self.setAbsoluteTolerance() 780 self.__adaptSubTolerance=adaptSubTolerance 781 self.verbose=False 782 def getSolverOptionsFlux(self): 783 """ 784 Returns the solver options used to solve the flux problems 785 786 *(I+D^*D)u=F* 787 788 :return: `SolverOptions` 789 """ 790 return self.__pde_v.getSolverOptions() 791 def setSolverOptionsFlux(self, options=None): 792 """ 793 Sets the solver options used to solve the flux problems 794 795 *(I+D^*D)u=F* 796 797 If ``options`` is not present, the options are reset to default 798 :param options: `SolverOptions` 799 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 800 """ 801 return self.__pde_v.setSolverOptions(options) 802 def getSolverOptionsPressure(self): 803 """ 804 Returns the solver options used to solve the pressure problems 805 806 *(Q^*Q)p=Q^*G* 807 808 :return: `SolverOptions` 809 """ 810 return self.__pde_p.getSolverOptions() 811 def setSolverOptionsPressure(self, options=None): 812 """ 813 Sets the solver options used to solve the pressure problems 814 815 *(Q^*Q)p=Q^*G* 816 817 If ``options`` is not present, the options are reset to default 818 :param options: `SolverOptions` 819 :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called. 820 """ 821 return self.__pde_p.setSolverOptions(options) 822 823 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): 824 """ 825 assigns values to model parameters 826 827 :param f: volumetic sources/sinks 828 :type f: scalar value on the domain (e.g. `escript.Data`) 829 :param g: flux sources/sinks 830 :type g: vector values on the domain (e.g. `escript.Data`) 831 :param location_of_fixed_pressure: mask for locations where pressure is fixed 832 :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`) 833 :param location_of_fixed_flux: mask for locations where flux is fixed. 834 :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`) 835 :param permeability: permeability tensor. If scalar ``s`` is given the tensor with 836 ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with 837 ``v`` on the main diagonal is used. 838 :type permeability: scalar, vector or tensor values on the domain (e.g. `escript.Data`) 839 840 :note: the values of parameters which are not set by calling ``setValue`` are not altered. 841 :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0) 842 or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal 843 is along the *x_i* axis. 844 """ 845 if f !=None: 846 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X")) 847 if f.isEmpty(): 848 f=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) 849 else: 850 if f.getRank()>0: raise ValueError,"illegal rank of f." 851 self.__f=f 852 if g !=None: 853 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y")) 854 if g.isEmpty(): 855 g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) 856 else: 857 if not g.getShape()==(self.domain.getDim(),): 858 raise ValueError,"illegal shape of g" 859 self.__g=g 860 861 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure) 862 if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux) 863 864 if permeability!=None: 865 perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A")) 866 if perm.getRank()==0: 867 perm=perm*util.kronecker(self.domain.getDim()) 868 elif perm.getRank()==1: 869 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm 870 for i in range(self.domain.getDim()): perm[i,i]=perm2[i] 871 elif perm.getRank()==2: 872 pass 873 else: 874 raise ValueError,"illegal rank of permeability." 875 self.__permeability=perm 876 self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)) 877 878 def setTolerance(self,rtol=1e-4): 879 """ 880 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if 881 882 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )* 883 884 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}*. 885 886 :param rtol: relative tolerance for the pressure 887 :type rtol: non-negative ``float`` 888 """ 889 if rtol<0: 890 raise ValueError,"Relative tolerance needs to be non-negative." 891 self.__rtol=rtol 892 def getTolerance(self): 893 """ 894 returns the relative tolerance 895 896 :return: current relative tolerance 897 :rtype: ``float`` 898 """ 899 return self.__rtol 900 901 def setAbsoluteTolerance(self,atol=0.): 902 """ 903 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if 904 905 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )* 906 907 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}*. 908 909 :param atol: absolute tolerance for the pressure 910 :type atol: non-negative ``float`` 911 """ 912 if atol<0: 913 raise ValueError,"Absolute tolerance needs to be non-negative." 914 self.__atol=atol 915 def getAbsoluteTolerance(self): 916 """ 917 returns the absolute tolerance 918 919 :return: current absolute tolerance 920 :rtype: ``float`` 921 """ 922 return self.__atol 923 def getSubProblemTolerance(self): 924 """ 925 Returns a suitable subtolerance 926 @type: ``float`` 927 """ 928 return max(util.EPSILON**(0.75),self.getTolerance()**2) 929 def setSubProblemTolerance(self): 930 """ 931 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected. 932 """ 933 if self.__adaptSubTolerance: 934 sub_tol=self.getSubProblemTolerance() 935 self.getSolverOptionsFlux().setTolerance(sub_tol) 936 self.getSolverOptionsFlux().setAbsoluteTolerance(0.) 937 self.getSolverOptionsPressure().setTolerance(sub_tol) 938 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 939 if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol 940 941 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10): 942 """ 943 solves the problem. 944 945 The iteration is terminated if the residual norm is less then self.getTolerance(). 946 947 :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. 948 :type u0: vector value on the domain (e.g. `escript.Data`). 949 :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. 950 :type p0: scalar value on the domain (e.g. `escript.Data`). 951 :param verbose: if set some information on iteration progress are printed 952 :type verbose: ``bool`` 953 :return: flux and pressure 954 :rtype: ``tuple`` of `escript.Data`. 955 956 :note: The problem is solved as a least squares form 957 958 *(I+D^*D)u+Qp=D^*f+g* 959 *Q^*u+Q^*Qp=Q^*g* 960 961 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*. 962 We eliminate the flux form the problem by setting 963 964 *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux 965 966 form the first equation. Inserted into the second equation we get 967 968 *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0 on location_of_fixed_pressure 969 970 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step 971 PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme. 972 """ 973 self.verbose=verbose 974 rtol=self.getTolerance() 975 atol=self.getAbsoluteTolerance() 976 self.setSubProblemTolerance() 977 num_corrections=0 978 converged=False 979 p=p0 980 norm_r=None 981 while not converged: 982 v=self.getFlux(p, fixed_flux=u0) 983 Qp=self.__Q(p) 984 norm_v=self.__L2(v) 985 norm_Qp=self.__L2(Qp) 986 if norm_v == 0.: 987 if norm_Qp == 0.: 988 return v,p 989 else: 990 fac=norm_Qp 991 else: 992 if norm_Qp == 0.: 993 fac=norm_v 994 else: 995 fac=2./(1./norm_v+1./norm_Qp) 996 ATOL=(atol+rtol*fac) 997 if self.verbose: 998 print "DarcyFlux: L2 norm of v = %e."%norm_v 999 print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp 1000 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,escript.Function(self.domain))-Qp)**2)**(0.5),) 1001 print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),) 1002 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL 1003 if norm_r == None or norm_r>ATOL: 1004 if num_corrections>max_num_corrections: 1005 raise ValueError,"maximum number of correction steps reached." 1006 p,r, norm_r=PCG(self.__g-util.interpolate(v,escript.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) 1007 num_corrections+=1 1008 else: 1009 converged=True 1010 return v,p 1011 def __L2(self,v): 1012 return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2)) 1013 1014 def __Q(self,p): 1015 return util.tensor_mult(self.__permeability,util.grad(p)) 1016 1017 def __Aprod(self,dp): 1018 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator" 1019 Qdp=self.__Q(dp) 1020 self.__pde_v.setValue(Y=-Qdp,X=escript.Data(), r=escript.Data()) 1021 du=self.__pde_v.getSolution() 1022 return Qdp+du 1023 def __inner_GMRES(self,r,s): 1024 return util.integrate(util.inner(r,s)) 1025 1026 def __inner_PCG(self,p,r): 1027 return util.integrate(util.inner(self.__Q(p), r)) 1028 1029 def __Msolve_PCG(self,r): 1030 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner" 1031 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=escript.Data(), r=escript.Data()) 1032 return self.__pde_p.getSolution() 1033 1034 def getFlux(self,p=None, fixed_flux=escript.Data()): 1035 """ 1036 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux`` 1037 on locations where ``location_of_fixed_flux`` is positive (see `setValue`). 1038 Note that ``g`` and ``f`` are used, see `setValue`. 1039 1040 :param p: pressure. 1041 :type p: scalar value on the domain (e.g. `escript.Data`). 1042 :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``. 1043 :type fixed_flux: vector values on the domain (e.g. `escript.Data`). 1044 :return: flux 1045 :rtype: `escript.Data` 1046 :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}* 1047 for the permeability *k_{ij}* 1048 """ 1049 self.setSubProblemTolerance() 1050 g=self.__g 1051 f=self.__f 1052 self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux) 1053 if p == None: 1054 self.__pde_v.setValue(Y=g) 1055 else: 1056 self.__pde_v.setValue(Y=g-self.__Q(p)) 1057 return self.__pde_v.getSolution() 1058 1059 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 1060 """ 1061 solves 1062 1063 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j} 1064 u_{i,i}=0 1065 1066 u=0 where fixed_u_mask>0 1067 eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j 1068 1069 if surface_stress is not given 0 is assumed. 1070 1071 typical usage: 1072 1073 sp=StokesProblemCartesian(domain) 1074 sp.setTolerance() 1075 sp.initialize(...) 1076 v,p=sp.solve(v0,p0) 1077 """ 1078 def __init__(self,domain,**kwargs): 1079 """ 1080 initialize the Stokes Problem 1081 1082 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be 1083 LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation 1084 with macro elements for the pressure. 1085 1086 :param domain: domain of the problem. 1087 :type domain: `Domain` 1088 """ 1089 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 1090 self.domain=domain 1091 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 1092 self.__pde_u.setSymmetryOn() 1093 1094 self.__pde_prec=LinearPDE(domain) 1095 self.__pde_prec.setReducedOrderOn() 1096 self.__pde_prec.setSymmetryOn() 1097 1098 self.__pde_proj=LinearPDE(domain) 1099 self.__pde_proj.setReducedOrderOn() 1100 self.__pde_proj.setValue(D=1) 1101 self.__pde_proj.setSymmetryOn() 1102 1103 def getSolverOptionsVelocity(self): 1104 """ 1105 returns the solver options used solve the equation for velocity. 1106 1107 :rtype: `SolverOptions` 1108 """ 1109 return self.__pde_u.getSolverOptions() 1110 def setSolverOptionsVelocity(self, options=None): 1111 """ 1112 set the solver options for solving the equation for velocity. 1113 1114 :param options: new solver options 1115 :type options: `SolverOptions` 1116 """ 1117 self.__pde_u.setSolverOptions(options) 1118 def getSolverOptionsPressure(self): 1119 """ 1120 returns the solver options used solve the equation for pressure. 1121 :rtype: `SolverOptions` 1122 """ 1123 return self.__pde_prec.getSolverOptions() 1124 def setSolverOptionsPressure(self, options=None): 1125 """ 1126 set the solver options for solving the equation for pressure. 1127 :param options: new solver options 1128 :type options: `SolverOptions` 1129 """ 1130 self.__pde_prec.setSolverOptions(options) 1131 1132 def setSolverOptionsDiv(self, options=None): 1133 """ 1134 set the solver options for solving the equation to project the divergence of 1135 the velocity onto the function space of presure. 1136 1137 :param options: new solver options 1138 :type options: `SolverOptions` 1139 """ 1140 self.__pde_proj.setSolverOptions(options) 1141 def getSolverOptionsDiv(self): 1142 """ 1143 returns the solver options for solving the equation to project the divergence of 1144 the velocity onto the function space of presure. 1145 1146 :rtype: `SolverOptions` 1147 """ 1148 return self.__pde_proj.getSolverOptions() 1149 1150 def updateStokesEquation(self, v, p): 1151 """ 1152 updates the Stokes equation to consider dependencies from ``v`` and ``p`` 1153 :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values. 1154 """ 1155 pass 1156 def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None): 1157 """ 1158 assigns new values to the model parameters. 1159 1160 :param f: external force 1161 :type f: `Vector` object in `FunctionSpace` `Function` or similar 1162 :param fixed_u_mask: mask of locations with fixed velocity. 1163 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar 1164 :param eta: viscosity 1165 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar 1166 :param surface_stress: normal surface stress 1167 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar 1168 :param stress: initial stress 1169 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar 1170 """ 1171 if eta !=None: 1172 k=util.kronecker(self.domain.getDim()) 1173 kk=util.outer(k,k) 1174 self.eta=util.interpolate(eta, escript.Function(self.domain)) 1175 self.__pde_prec.setValue(D=1/self.eta) 1176 self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3))) 1177 if restoration_factor!=None: 1178 n=self.domain.getNormal() 1179 self.__pde_u.setValue(d=restoration_factor*util.outer(n,n)) 1180 if fixed_u_mask!=None: 1181 self.__pde_u.setValue(q=fixed_u_mask) 1182 if f!=None: self.__f=f 1183 if surface_stress!=None: self.__surface_stress=surface_stress 1184 if stress!=None: self.__stress=stress 1185 1186 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0): 1187 """ 1188 assigns values to the model parameters 1189 1190 :param f: external force 1191 :type f: `Vector` object in `FunctionSpace` `Function` or similar 1192 :param fixed_u_mask: mask of locations with fixed velocity. 1193 :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar 1194 :param eta: viscosity 1195 :type eta: `Scalar` object on `FunctionSpace` `Function` or similar 1196 :param surface_stress: normal surface stress 1197 :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar 1198 :param stress: initial stress 1199 :type stress: `Tensor` object on `FunctionSpace` `Function` or similar 1200 """ 1201 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor) 1202 1203 def Bv(self,v,tol): 1204 """ 1205 returns inner product of element p and div(v) 1206 1207 :param v: a residual 1208 :return: inner product of element p and div(v) 1209 :rtype: ``float`` 1210 """ 1211 self.__pde_proj.setValue(Y=-util.div(v)) 1212 self.getSolverOptionsDiv().setTolerance(tol) 1213 self.getSolverOptionsDiv().setAbsoluteTolerance(0.) 1214 out=self.__pde_proj.getSolution() 1215 return out 1216 1217 def inner_pBv(self,p,Bv): 1218 """ 1219 returns inner product of element p and Bv=-div(v) 1220 1221 :param p: a pressure increment 1222 :param Bv: a residual 1223 :return: inner product of element p and Bv=-div(v) 1224 :rtype: ``float`` 1225 """ 1226 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain))) 1227 1228 def inner_p(self,p0,p1): 1229 """ 1230 Returns inner product of p0 and p1 1231 1232 :param p0: a pressure 1233 :param p1: a pressure 1234 :return: inner product of p0 and p1 1235 :rtype: ``float`` 1236 """ 1237 s0=util.interpolate(p0, escript.Function(self.domain)) 1238 s1=util.interpolate(p1, escript.Function(self.domain)) 1239 return util.integrate(s0*s1) 1240 1241 def norm_v(self,v): 1242 """ 1243 returns the norm of v 1244 1245 :param v: a velovity 1246 :return: norm of v 1247 :rtype: non-negative ``float`` 1248 """ 1249 return util.sqrt(util.integrate(util.length(util.grad(v))**2)) 1250 1251 1252 def getDV(self, p, v, tol): 1253 """ 1254 return the value for v for a given p (overwrite) 1255 1256 :param p: a pressure 1257 :param v: a initial guess for the value v to return. 1258 :return: dv given as *Adv=(f-Av-B^*p)* 1259 """ 1260 self.updateStokesEquation(v,p) 1261 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress) 1262 self.getSolverOptionsVelocity().setTolerance(tol) 1263 self.getSolverOptionsVelocity().setAbsoluteTolerance(0.) 1264 if self.__stress.isEmpty(): 1265 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 1266 else: 1267 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v))) 1268 out=self.__pde_u.getSolution() 1269 return out 1270 1271 def norm_Bv(self,Bv): 1272 """ 1273 Returns Bv (overwrite). 1274 1275 :rtype: equal to the type of p 1276 :note: boundary conditions on p should be zero! 1277 """ 1278 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2)) 1279 1280 def solve_AinvBt(self,p, tol): 1281 """ 1282 Solves *Av=B^*p* with accuracy `tol` 1283 1284 :param p: a pressure increment 1285 :return: the solution of *Av=B^*p* 1286 :note: boundary conditions on v should be zero! 1287 """ 1288 self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain)) 1289 out=self.__pde_u.getSolution() 1290 return out 1291 1292 def solve_prec(self,Bv, tol): 1293 """ 1294 applies preconditioner for for *BA^{-1}B^** to *Bv* 1295 with accuracy `self.getSubProblemTolerance()` 1296 1297 :param Bv: velocity increment 1298 :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )* 1299 :note: boundary conditions on p are zero. 1300 """ 1301 self.__pde_prec.setValue(Y=Bv) 1302 self.getSolverOptionsPressure().setTolerance(tol) 1303 self.getSolverOptionsPressure().setAbsoluteTolerance(0.) 1304 out=self.__pde_prec.getSolution() 1305 return out

 ViewVC Help Powered by ViewVC 1.1.26