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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26