Contents of /branches/scons-dev/escript/py_src/flows.py

Revision 1681 - (show annotations)
Thu Jul 31 06:20:18 2008 UTC (11 years, 8 months ago) by ksteube
File MIME type: text/x-python
File size: 9782 byte(s)
```Cleaned up some problems from the automatic merge

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

 ViewVC Help Powered by ViewVC 1.1.26