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

Revision 1465 - (show annotations)
Wed Apr 2 03:28:25 2008 UTC (11 years, 4 months ago) by artak
Original Path: trunk/escript/py_src/flows.py
File MIME type: text/x-python
File size: 4697 byte(s)

 1 # \$Id:\$ 2 # 3 ####################################################### 4 # 5 # Copyright 2008 by University of Queensland 6 # 7 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ####################################################### 13 # 14 15 """ 16 Some models for flow 17 18 @var __author__: name of author 19 @var __copyright__: copyrights 20 @var __license__: licence agreement 21 @var __url__: url entry point on documentation 22 @var __version__: version 23 @var __date__: date of the version 24 """ 25 26 __author__="Lutz Gross, l.gross@uq.edu.au" 27 __copyright__=""" Copyright (c) 2008 by ACcESS MNRF 28 29 Primary Business: Queensland, Australia""" 30 __license__="""Licensed under the Open Software License version 3.0 31 32 __url__= 33 __version__="\$Revision:\$" 34 __date__="\$Date:\$" 35 36 from escript import * 37 import util 38 from linearPDEs import LinearPDE 39 from pdetools import HomogeneousSaddlePointProblem 40 41 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 42 """ 43 solves 44 45 -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i 46 u_{i,i}=0 47 48 u=0 where fixed_u_mask>0 49 eta*(u_{i,j}+u_{j,i})*n_j=surface_stress 50 51 if surface_stress is not give 0 is assumed. 52 53 typical usage: 54 55 sp=StokesProblemCartesian(domain) 56 sp.setTolerance() 57 sp.initialize(...) 58 v,p=sp.solve(v0,p0) 59 """ 60 def __init__(self,domain,**kwargs): 61 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 62 self.domain=domain 63 self.vol=util.integrate(1.,Function(self.domain)) 64 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 65 self.__pde_u.setSymmetryOn() 66 self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0) 67 68 self.__pde_prec=LinearPDE(domain) 69 self.__pde_prec.setReducedOrderOn() 70 self.__pde_prec.setSymmetryOn() 71 72 self.__pde_proj=LinearPDE(domain) 73 self.__pde_proj.setReducedOrderOn() 74 self.__pde_proj.setSymmetryOn() 75 self.__pde_proj.setValue(D=1.) 76 77 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()): 78 self.eta=eta 79 A =self.__pde_u.createCoefficientOfGeneralPDE("A") 80 self.__pde_u.setValue(A=Data()) 81 for i in range(self.domain.getDim()): 82 for j in range(self.domain.getDim()): 83 A[i,j,j,i] += 1. 84 A[i,j,i,j] += 1. 85 self.__pde_prec.setValue(D=1./self.eta) 86 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) 87 88 def B(self,arg): 89 d=util.div(arg) 90 self.__pde_proj.setValue(Y=d) 91 self.__pde_proj.setTolerance(self.getSubProblemTolerance()) 92 return self.__pde_proj.getSolution(verbose=self.show_details) 93 94 def inner(self,p0,p1): 95 s0=util.interpolate(p0,Function(self.domain)) 96 s1=util.interpolate(p1,Function(self.domain)) 97 return util.integrate(s0*s1) 98 99 def getStress(self,u): 100 mg=util.grad(u) 101 return 2.*self.eta*util.symmetric(mg) 102 103 def solve_A(self,u,p): 104 """ 105 solves Av=f-Au-B^*p (v=0 on fixed_u_mask) 106 """ 107 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 108 self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain)) 109 return self.__pde_u.getSolution(verbose=self.show_details) 110 111 def solve_prec(self,p): 112 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 113 self.__pde_prec.setValue(Y=p) 114 q=self.__pde_prec.getSolution(verbose=self.show_details) 115 return q 116 def stoppingcriterium(self,Bv,v,p): 117 n_r=util.sqrt(self.inner(Bv,Bv)) 118 n_v=util.Lsup(v) 119 if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v) 120 self.iter+=1 121 if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance(): 122 if self.verbose: print "PCG terminated after %s steps."%self.iter 123 return True 124 else: 125 return False 126 def stoppingcriterium_GMRES(self,rho,r): 127 if self.verbose: print "GMRES step %s: L2(rho) = %s, L2(b)*TOL=%s"%(self.iter,rho,r*self.getTolerance()) 128 self.iter+=1 129 if rho <= r*self.getTolerance(): 130 if self.verbose: print "GMRES terminated after %s steps."%self.iter 131 return True 132 else: 133 return False

 ViewVC Help Powered by ViewVC 1.1.26