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

  ViewVC Help
Powered by ViewVC 1.1.26