/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3452 - (show annotations)
Tue Jan 25 01:53:57 2011 UTC (8 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 44103 byte(s)
Fixed building cookbook, including toc (scons scans for keywords in the tex
files but not in .cls as it seems).
Also fixed a few epydoc warnings and typos.


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

  ViewVC Help
Powered by ViewVC 1.1.26