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

Revision 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years, 3 months ago) by ksteube
File MIME type: text/x-python
File size: 9794 byte(s)
```Copyright updated in all python files

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2008 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-2008 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 from escript import * 36 import util 37 from linearPDEs import LinearPDE 38 from pdetools import HomogeneousSaddlePointProblem,Projector 39 40 class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem): 41 """ 42 solves 43 44 -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i 45 u_{i,i}=0 46 47 u=0 where fixed_u_mask>0 48 eta*(u_{i,j}+u_{j,i})*n_j=surface_stress 49 50 if surface_stress is not give 0 is assumed. 51 52 typical usage: 53 54 sp=StokesProblemCartesian(domain) 55 sp.setTolerance() 56 sp.initialize(...) 57 v,p=sp.solve(v0,p0) 58 """ 59 def __init__(self,domain,**kwargs): 60 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 61 self.domain=domain 62 self.vol=util.integrate(1.,Function(self.domain)) 63 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 64 self.__pde_u.setSymmetryOn() 65 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0) 66 67 # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1) 68 # self.__pde_proj.setReducedOrderOn() 69 # self.__pde_proj.setSymmetryOn() 70 # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING) 71 72 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()): 73 self.eta=eta 74 A =self.__pde_u.createCoefficientOfGeneralPDE("A") 75 self.__pde_u.setValue(A=Data()) 76 for i in range(self.domain.getDim()): 77 for j in range(self.domain.getDim()): 78 A[i,j,j,i] += 1. 79 A[i,j,i,j] += 1. 80 # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain)) 81 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) 82 83 # self.__pde_proj.setValue(D=1/eta) 84 # self.__pde_proj.setValue(Y=1.) 85 # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain)) 86 self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain)) 87 88 def B(self,arg): 89 a=util.div(arg, ReducedFunction(self.domain)) 90 return a-util.integrate(a)/self.vol 91 92 def inner(self,p0,p1): 93 return util.integrate(p0*p1) 94 # return util.integrate(1/self.__inv_eta**2*p0*p1) 95 96 def getStress(self,u): 97 mg=util.grad(u) 98 return 2.*self.eta*util.symmetric(mg) 99 def getEtaEffective(self): 100 return self.eta 101 102 def solve_A(self,u,p): 103 """ 104 solves Av=f-Au-B^*p (v=0 on fixed_u_mask) 105 """ 106 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 107 self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain)) 108 return self.__pde_u.getSolution(verbose=self.show_details) 109 110 111 def solve_prec(self,p): 112 a=self.__inv_eta*p 113 return a-util.integrate(a)/self.vol 114 115 def stoppingcriterium(self,Bv,v,p): 116 n_r=util.sqrt(self.inner(Bv,Bv)) 117 n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2)) 118 if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v) 119 if self.iter == 0: self.__n_v=n_v; 120 self.__n_v, n_v_old =n_v, self.__n_v 121 self.iter+=1 122 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance(): 123 if self.verbose: print "PCG terminated after %s steps."%self.iter 124 return True 125 else: 126 return False 127 def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None): 128 if TOL==None: 129 TOL=self.getTolerance() 130 if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL) 131 self.iter+=1 132 133 if norm_r <= norm_b*TOL: 134 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter) 135 return True 136 else: 137 return False 138 139 140 class StokesProblemCartesian(HomogeneousSaddlePointProblem): 141 """ 142 solves 143 144 -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i 145 u_{i,i}=0 146 147 u=0 where fixed_u_mask>0 148 eta*(u_{i,j}+u_{j,i})*n_j=surface_stress 149 150 if surface_stress is not give 0 is assumed. 151 152 typical usage: 153 154 sp=StokesProblemCartesian(domain) 155 sp.setTolerance() 156 sp.initialize(...) 157 v,p=sp.solve(v0,p0) 158 """ 159 def __init__(self,domain,**kwargs): 160 HomogeneousSaddlePointProblem.__init__(self,**kwargs) 161 self.domain=domain 162 self.vol=util.integrate(1.,Function(self.domain)) 163 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) 164 self.__pde_u.setSymmetryOn() 165 # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0) 166 167 self.__pde_prec=LinearPDE(domain) 168 self.__pde_prec.setReducedOrderOn() 169 self.__pde_prec.setSymmetryOn() 170 171 self.__pde_proj=LinearPDE(domain) 172 self.__pde_proj.setReducedOrderOn() 173 self.__pde_proj.setSymmetryOn() 174 self.__pde_proj.setValue(D=1.) 175 176 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()): 177 self.eta=eta 178 A =self.__pde_u.createCoefficientOfGeneralPDE("A") 179 self.__pde_u.setValue(A=Data()) 180 for i in range(self.domain.getDim()): 181 for j in range(self.domain.getDim()): 182 A[i,j,j,i] += 1. 183 A[i,j,i,j] += 1. 184 self.__pde_prec.setValue(D=1/self.eta) 185 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) 186 187 def B(self,arg): 188 d=util.div(arg) 189 self.__pde_proj.setValue(Y=d) 190 self.__pde_proj.setTolerance(self.getSubProblemTolerance()) 191 return self.__pde_proj.getSolution(verbose=self.show_details) 192 193 def inner(self,p0,p1): 194 s0=util.interpolate(p0,Function(self.domain)) 195 s1=util.interpolate(p1,Function(self.domain)) 196 return util.integrate(s0*s1) 197 198 def inner_a(self,a0,a1): 199 p0=util.interpolate(a0[1],Function(self.domain)) 200 p1=util.interpolate(a1[1],Function(self.domain)) 201 alfa=(1/self.vol)*util.integrate(p0) 202 beta=(1/self.vol)*util.integrate(p1) 203 v0=util.grad(a0[0]) 204 v1=util.grad(a1[0]) 205 return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1)) 206 207 208 def getStress(self,u): 209 mg=util.grad(u) 210 return 2.*self.eta*util.symmetric(mg) 211 def getEtaEffective(self): 212 return self.eta 213 214 def solve_A(self,u,p): 215 """ 216 solves Av=f-Au-B^*p (v=0 on fixed_u_mask) 217 """ 218 self.__pde_u.setTolerance(self.getSubProblemTolerance()) 219 self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain)) 220 return self.__pde_u.getSolution(verbose=self.show_details) 221 222 223 def solve_prec(self,p): 224 #proj=Projector(domain=self.domain, reduce = True, fast=False) 225 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 226 self.__pde_prec.setValue(Y=p) 227 q=self.__pde_prec.getSolution(verbose=self.show_details) 228 return q 229 230 def solve_prec1(self,p): 231 #proj=Projector(domain=self.domain, reduce = True, fast=False) 232 self.__pde_prec.setTolerance(self.getSubProblemTolerance()) 233 self.__pde_prec.setValue(Y=p) 234 q=self.__pde_prec.getSolution(verbose=self.show_details) 235 q0=util.interpolate(q,Function(self.domain)) 236 print util.inf(q*q0),util.sup(q*q0) 237 q-=(1/self.vol)*util.integrate(q0) 238 print util.inf(q*q0),util.sup(q*q0) 239 return q 240 241 def stoppingcriterium(self,Bv,v,p): 242 n_r=util.sqrt(self.inner(Bv,Bv)) 243 n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2)) 244 if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) 245 if self.iter == 0: self.__n_v=n_v; 246 self.__n_v, n_v_old =n_v, self.__n_v 247 self.iter+=1 248 if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance(): 249 if self.verbose: print "PCG terminated after %s steps."%self.iter 250 return True 251 else: 252 return False 253 def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None): 254 if TOL==None: 255 TOL=self.getTolerance() 256 if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL) 257 self.iter+=1 258 259 if norm_r <= norm_b*TOL: 260 if self.verbose: print "%s terminated after %s steps."%(solver,self.iter) 261 return True 262 else: 263 return False 264 265

 ViewVC Help Powered by ViewVC 1.1.26