1 |
# $Id:$ |
2 |
# |
3 |
####################################################### |
4 |
# |
5 |
# Copyright 2008 by University of Queensland |
6 |
# |
7 |
# http://esscc.uq.edu.au |
8 |
# Primary Business: Queensland, Australia |
9 |
# Licensed under the Open Software License version 3.0 |
10 |
# http://www.opensource.org/licenses/osl-3.0.php |
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 |
http://www.access.edu.au |
29 |
Primary Business: Queensland, Australia""" |
30 |
__license__="""Licensed under the Open Software License version 3.0 |
31 |
http://www.opensource.org/licenses/osl-3.0.php""" |
32 |
__url__="http://www.iservo.edu.au/esys" |
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 |
|
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 |
|
128 |
|
129 |
class StokesProblemCartesian(HomogeneousSaddlePointProblem): |
130 |
""" |
131 |
solves |
132 |
|
133 |
-(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i |
134 |
u_{i,i}=0 |
135 |
|
136 |
u=0 where fixed_u_mask>0 |
137 |
eta*(u_{i,j}+u_{j,i})*n_j=surface_stress |
138 |
|
139 |
if surface_stress is not give 0 is assumed. |
140 |
|
141 |
typical usage: |
142 |
|
143 |
sp=StokesProblemCartesian(domain) |
144 |
sp.setTolerance() |
145 |
sp.initialize(...) |
146 |
v,p=sp.solve(v0,p0) |
147 |
""" |
148 |
def __init__(self,domain,**kwargs): |
149 |
HomogeneousSaddlePointProblem.__init__(self,**kwargs) |
150 |
self.domain=domain |
151 |
self.vol=util.integrate(1.,Function(self.domain)) |
152 |
self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim()) |
153 |
self.__pde_u.setSymmetryOn() |
154 |
# self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0) |
155 |
|
156 |
self.__pde_prec=LinearPDE(domain) |
157 |
self.__pde_prec.setReducedOrderOn() |
158 |
self.__pde_prec.setSymmetryOn() |
159 |
|
160 |
self.__pde_proj=LinearPDE(domain) |
161 |
self.__pde_proj.setReducedOrderOn() |
162 |
self.__pde_proj.setSymmetryOn() |
163 |
self.__pde_proj.setValue(D=1.) |
164 |
|
165 |
def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()): |
166 |
self.eta=eta |
167 |
A =self.__pde_u.createCoefficientOfGeneralPDE("A") |
168 |
self.__pde_u.setValue(A=Data()) |
169 |
for i in range(self.domain.getDim()): |
170 |
for j in range(self.domain.getDim()): |
171 |
A[i,j,j,i] += 1. |
172 |
A[i,j,i,j] += 1. |
173 |
self.__pde_prec.setValue(D=1/self.eta) |
174 |
self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress) |
175 |
|
176 |
def B(self,arg): |
177 |
d=util.div(arg) |
178 |
self.__pde_proj.setValue(Y=d) |
179 |
self.__pde_proj.setTolerance(self.getSubProblemTolerance()) |
180 |
return self.__pde_proj.getSolution(verbose=self.show_details) |
181 |
|
182 |
def inner(self,p0,p1): |
183 |
s0=util.interpolate(p0,Function(self.domain)) |
184 |
s1=util.interpolate(p1,Function(self.domain)) |
185 |
return util.integrate(s0*s1) |
186 |
|
187 |
def inner_a(self,a0,a1): |
188 |
p0=util.interpolate(a0[1],Function(self.domain)) |
189 |
p1=util.interpolate(a1[1],Function(self.domain)) |
190 |
alfa=(1/self.vol)*util.integrate(p0) |
191 |
beta=(1/self.vol)*util.integrate(p1) |
192 |
v0=util.grad(a0[0]) |
193 |
v1=util.grad(a1[0]) |
194 |
return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1)) |
195 |
|
196 |
|
197 |
def getStress(self,u): |
198 |
mg=util.grad(u) |
199 |
return 2.*self.eta*util.symmetric(mg) |
200 |
def getEtaEffective(self): |
201 |
return self.eta |
202 |
|
203 |
def solve_A(self,u,p): |
204 |
""" |
205 |
solves Av=f-Au-B^*p (v=0 on fixed_u_mask) |
206 |
""" |
207 |
self.__pde_u.setTolerance(self.getSubProblemTolerance()) |
208 |
self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain)) |
209 |
return self.__pde_u.getSolution(verbose=self.show_details) |
210 |
|
211 |
|
212 |
def solve_prec(self,p): |
213 |
#proj=Projector(domain=self.domain, reduce = True, fast=False) |
214 |
self.__pde_prec.setTolerance(self.getSubProblemTolerance()) |
215 |
self.__pde_prec.setValue(Y=p) |
216 |
q=self.__pde_prec.getSolution(verbose=self.show_details) |
217 |
return q |
218 |
|
219 |
def solve_prec1(self,p): |
220 |
#proj=Projector(domain=self.domain, reduce = True, fast=False) |
221 |
self.__pde_prec.setTolerance(self.getSubProblemTolerance()) |
222 |
self.__pde_prec.setValue(Y=p) |
223 |
q=self.__pde_prec.getSolution(verbose=self.show_details) |
224 |
q0=util.interpolate(q,Function(self.domain)) |
225 |
print util.inf(q*q0),util.sup(q*q0) |
226 |
q-=(1/self.vol)*util.integrate(q0) |
227 |
print util.inf(q*q0),util.sup(q*q0) |
228 |
return q |
229 |
|
230 |
def stoppingcriterium(self,Bv,v,p): |
231 |
n_r=util.sqrt(self.inner(Bv,Bv)) |
232 |
n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2)) |
233 |
if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) |
234 |
if self.iter == 0: self.__n_v=n_v; |
235 |
self.__n_v, n_v_old =n_v, self.__n_v |
236 |
self.iter+=1 |
237 |
if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance(): |
238 |
if self.verbose: print "PCG terminated after %s steps."%self.iter |
239 |
return True |
240 |
else: |
241 |
return False |
242 |
def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None): |
243 |
if TOL==None: |
244 |
TOL=self.getTolerance() |
245 |
if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL) |
246 |
self.iter+=1 |
247 |
|
248 |
if norm_r <= norm_b*TOL: |
249 |
if self.verbose: print "%s terminated after %s steps."%(solver,self.iter) |
250 |
return True |
251 |
else: |
252 |
return False |
253 |
|
254 |
|