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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6414 - (show annotations)
Wed Nov 2 02:44:31 2016 UTC (2 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 69205 byte(s)
work around for bug #391
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2016 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import print_function, division
18
19 __copyright__="""Copyright (c) 2003-2016 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 """
27 Provides some tools related to PDEs.
28
29 Currently includes:
30 - Projector - to project a discontinuous function onto a continuous function
31 - Locator - to trace values in data objects at a certain location
32 - TimeIntegrationManager - to handle extrapolation in time
33 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
34
35 :var __author__: name of author
36 :var __copyright__: copyrights
37 :var __license__: licence agreement
38 :var __url__: url entry point on documentation
39 :var __version__: version
40 :var __date__: date of the version
41 """
42
43 __author__="Lutz Gross, l.gross@uq.edu.au"
44
45
46 from . import escriptcpp as escore
47 from . import linearPDEs
48 from . import util
49 import math
50 import numpy
51
52 class TimeIntegrationManager(object):
53 """
54 A simple mechanism to manage time dependend values.
55
56 Typical usage is::
57
58 dt=0.1 # time increment
59 tm=TimeIntegrationManager(inital_value,p=1)
60 while t<1.
61 v_guess=tm.extrapolate(dt) # extrapolate to t+dt
62 v=...
63 tm.checkin(dt,v)
64 t+=dt
65
66 :note: currently only p=1 is supported.
67 """
68 def __init__(self,*inital_values,**kwargs):
69 """
70 Sets up the value manager where ``inital_values`` are the initial values
71 and p is the order used for extrapolation.
72 """
73 if "p" in kwargs:
74 self.__p=kwargs["p"]
75 else:
76 self.__p=1
77 if "time" in kwargs:
78 self.__t=kwargs["time"]
79 else:
80 self.__t=0.
81 self.__v_mem=[inital_values]
82 self.__order=0
83 self.__dt_mem=[]
84 self.__num_val=len(inital_values)
85
86 def getTime(self):
87 return self.__t
88
89 def getValue(self):
90 out=self.__v_mem[0]
91 if len(out)==1:
92 return out[0]
93 else:
94 return out
95
96 def checkin(self,dt,*values):
97 """
98 Adds new values to the manager. The p+1 last values are lost.
99 """
100 o=min(self.__order+1,self.__p)
101 self.__order=min(self.__order+1,self.__p)
102 v_mem_new=[values]
103 dt_mem_new=[dt]
104 for i in range(o-1):
105 v_mem_new.append(self.__v_mem[i])
106 dt_mem_new.append(self.__dt_mem[i])
107 v_mem_new.append(self.__v_mem[o-1])
108 self.__order=o
109 self.__v_mem=v_mem_new
110 self.__dt_mem=dt_mem_new
111 self.__t+=dt
112
113 def extrapolate(self,dt):
114 """
115 Extrapolates to ``dt`` forward in time.
116 """
117 if self.__order==0:
118 out=self.__v_mem[0]
119 else:
120 out=[]
121 for i in range(self.__num_val):
122 out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
123
124 if len(out)==0:
125 return None
126 elif len(out)==1:
127 return out[0]
128 else:
129 return out
130
131
132 class Projector(object):
133 """
134 The Projector is a factory which projects a discontinuous function onto a
135 continuous function on a given domain.
136 """
137 def __init__(self, domain, reduce=True, fast=True):
138 """
139 Creates a continuous function space projector for a domain.
140
141 :param domain: Domain of the projection.
142 :param reduce: Flag to reduce projection order
143 :param fast: Flag to use a fast method based on matrix lumping
144 """
145 self.__pde = linearPDEs.LinearPDE(domain)
146 if fast:
147 self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
148 self.__pde.setSymmetryOn()
149 self.__pde.setReducedOrderTo(reduce)
150 self.__pde.setValue(D = 1.)
151 return
152 def getSolverOptions(self):
153 """
154 Returns the solver options of the PDE solver.
155
156 :rtype: `linearPDEs.SolverOptions`
157 """
158 return self.__pde.getSolverOptions()
159
160 def getValue(self, input_data):
161 """
162 Projects ``input_data`` onto a continuous function.
163
164 :param input_data: the data to be projected
165 """
166 return self(input_data)
167
168 def __call__(self, input_data):
169 """
170 Projects ``input_data`` onto a continuous function.
171
172 :param input_data: the data to be projected
173 """
174 out=escore.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
175 self.__pde.setValue(Y = escore.Data(), Y_reduced = escore.Data())
176 if input_data.getRank()==0:
177 self.__pde.setValue(Y = input_data)
178 out=self.__pde.getSolution()
179 elif input_data.getRank()==1:
180 for i0 in range(input_data.getShape()[0]):
181 self.__pde.setValue(Y = input_data[i0])
182 out[i0]=self.__pde.getSolution()
183 elif input_data.getRank()==2:
184 for i0 in range(input_data.getShape()[0]):
185 for i1 in range(input_data.getShape()[1]):
186 self.__pde.setValue(Y = input_data[i0,i1])
187 out[i0,i1]=self.__pde.getSolution()
188 elif input_data.getRank()==3:
189 for i0 in range(input_data.getShape()[0]):
190 for i1 in range(input_data.getShape()[1]):
191 for i2 in range(input_data.getShape()[2]):
192 self.__pde.setValue(Y = input_data[i0,i1,i2])
193 out[i0,i1,i2]=self.__pde.getSolution()
194 else:
195 for i0 in range(input_data.getShape()[0]):
196 for i1 in range(input_data.getShape()[1]):
197 for i2 in range(input_data.getShape()[2]):
198 for i3 in range(input_data.getShape()[3]):
199 self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
200 out[i0,i1,i2,i3]=self.__pde.getSolution()
201 return out
202
203 class NoPDE(object):
204 """
205 Solves the following problem for u:
206
207 *kronecker[i,j]*D[j]*u[j]=Y[i]*
208
209 with constraint
210
211 *u[j]=r[j]* where *q[j]>0*
212
213 where *D*, *Y*, *r* and *q* are given functions of rank 1.
214
215 In the case of scalars this takes the form
216
217 *D*u=Y*
218
219 with constraint
220
221 *u=r* where *q>0*
222
223 where *D*, *Y*, *r* and *q* are given scalar functions.
224
225 The constraint overwrites any other condition.
226
227 :note: This class is similar to the `linearPDEs.LinearPDE` class with
228 A=B=C=X=0 but has the intention that all input parameters are given
229 in `Solution` or `ReducedSolution`.
230 """
231 # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
232 # this.
233 def __init__(self,domain,D=None,Y=None,q=None,r=None):
234 """
235 Initializes the problem.
236
237 :param domain: domain of the PDE
238 :type domain: `Domain`
239 :param D: coefficient of the solution
240 :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
241 :param Y: right hand side
242 :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
243 :param q: location of constraints
244 :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
245 :param r: value of solution at locations of constraints
246 :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
247 """
248 self.__domain=domain
249 self.__D=D
250 self.__Y=Y
251 self.__q=q
252 self.__r=r
253 self.__u=None
254 self.__function_space=escore.Solution(self.__domain)
255
256 def setReducedOn(self):
257 """
258 Sets the `FunctionSpace` of the solution to `ReducedSolution`.
259 """
260 self.__function_space=escore.ReducedSolution(self.__domain)
261 self.__u=None
262
263 def setReducedOff(self):
264 """
265 Sets the `FunctionSpace` of the solution to `Solution`.
266 """
267 self.__function_space=escore.Solution(self.__domain)
268 self.__u=None
269
270 def setValue(self,D=None,Y=None,q=None,r=None):
271 """
272 Assigns values to the parameters.
273
274 :param D: coefficient of the solution
275 :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
276 :param Y: right hand side
277 :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
278 :param q: location of constraints
279 :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
280 :param r: value of solution at locations of constraints
281 :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
282 """
283 if not D is None:
284 self.__D=D
285 self.__u=None
286 if not Y is None:
287 self.__Y=Y
288 self.__u=None
289 if not q is None:
290 self.__q=q
291 self.__u=None
292 if not r is None:
293 self.__r=r
294 self.__u=None
295
296 def getSolution(self):
297 """
298 Returns the solution.
299
300 :return: the solution of the problem
301 :rtype: `Data` object in the `FunctionSpace` `Solution` or
302 `ReducedSolution`
303 """
304 if self.__u is None:
305 if self.__D is None:
306 raise ValueError("coefficient D is undefined")
307 D=escore.Data(self.__D,self.__function_space)
308 if D.getRank()>1:
309 raise ValueError("coefficient D must have rank 0 or 1")
310 if self.__Y is None:
311 self.__u=escore.Data(0.,D.getShape(),self.__function_space)
312 else:
313 self.__u=1./D*self.__Y
314 if not self.__q is None:
315 q=util.wherePositive(escore.Data(self.__q,self.__function_space))
316 self.__u*=(1.-q)
317 if not self.__r is None: self.__u+=q*self.__r
318 return self.__u
319
320 class Locator(object):
321 """
322 Locator provides access to the values of data objects at a given spatial
323 coordinate x.
324
325 In fact, a Locator object finds the sample in the set of samples of a
326 given function space or domain which is closest to the given point x.
327 """
328
329 def __init__(self,where,x=numpy.zeros((3,))):
330 """
331 Initializes a Locator to access values in Data objects on the Doamin
332 or FunctionSpace for the sample point which is closest to the given
333 point x.
334
335 :param where: function space
336 :type where: `escript.FunctionSpace`
337 :param x: location(s) of the Locator
338 :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
339 """
340 if isinstance(where,escore.FunctionSpace):
341 self.__function_space=where
342 else:
343 self.__function_space=escore.ContinuousFunction(where)
344 iterative=False
345 if isinstance(x, list):
346 if len(x)==0:
347 raise ValueError("At least one point must be given.")
348 try:
349 iter(x[0])
350 iterative=True
351 except TypeError:
352 iterative=False
353 xxx=self.__function_space.getX()
354 if iterative:
355 self.__id=[]
356 for p in x:
357 self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).internal_minGlobalDataPoint())
358 else:
359 self.__id=util.length(xxx-x[:self.__function_space.getDim()]).internal_minGlobalDataPoint()
360
361 def __str__(self):
362 """
363 Returns the coordinates of the Locator as a string.
364 """
365 x=self.getX()
366 if isinstance(x,list):
367 out="["
368 first=True
369 for xx in x:
370 if not first:
371 out+=","
372 else:
373 first=False
374 out+=str(xx)
375 out+="]>"
376 else:
377 out=str(x)
378 return out
379
380 def getX(self):
381 """
382 Returns the exact coordinates of the Locator.
383 """
384 return self(self.getFunctionSpace().getX())
385
386 def getFunctionSpace(self):
387 """
388 Returns the function space of the Locator.
389 """
390 return self.__function_space
391
392 def getId(self,item=None):
393 """
394 Returns the identifier of the location.
395 """
396 if item is None:
397 return self.__id
398 else:
399 if isinstance(self.__id,list):
400 return self.__id[item]
401 else:
402 return self.__id
403
404
405 def __call__(self,data):
406 """
407 Returns the value of data at the Locator of a Data object.
408 """
409 return self.getValue(data)
410
411 def getValue(self,data):
412 """
413 Returns the value of ``data`` at the Locator if ``data`` is a `Data`
414 object otherwise the object is returned.
415 """
416 if isinstance(data,escore.Data):
417 dat=util.interpolate(data,self.getFunctionSpace())
418 ii=self.getId()
419 r=data.getRank()
420 if isinstance(ii,list):
421 out=[]
422 for i in ii:
423 # workaround for bug #391
424 if dat.isComplex():
425 o=numpy.array(dat.real().getTupleForGlobalDataPoint(*i))+1j*numpy.array(dat.imag().getTupleForGlobalDataPoint(*i))
426 else:
427 o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
428 if data.getRank()==0:
429 out.append(o[0])
430 else:
431 out.append(o)
432 return out
433 else:
434 # workaround for bug #391
435 if dat.isComplex():
436 print "DDD"
437 out=numpy.array(dat.real().getTupleForGlobalDataPoint(*ii))+1j*numpy.array(dat.imag().getTupleForGlobalDataPoint(*ii))
438 else:
439 out=numpy.array(dat.getTupleForGlobalDataPoint(*ii))
440 if data.getRank()==0:
441 return out[0]
442 else:
443 return out
444 else:
445 return data
446
447 def setValue(self, data, v):
448 """
449 Sets the value of the ``data`` at the Locator.
450 """
451 if isinstance(data, escore.Data):
452 if data.getFunctionSpace()!=self.getFunctionSpace():
453 raise TypeError("setValue: FunctionSpace of Locator and Data object must match.")
454 data.expand()
455 id=self.getId()
456 if isinstance(id, list):
457 for i in id:
458 data._setTupleForGlobalDataPoint(i[1], i[0], v)
459 else:
460 data._setTupleForGlobalDataPoint(id[1], id[0], v)
461 else:
462 raise TypeError("setValue: Invalid argument type.")
463
464
465 def getInfLocator(arg):
466 """
467 Return a Locator for a point with the inf value over all arg.
468 """
469 if not isinstance(arg, escore.Data):
470 raise TypeError("getInfLocator: Unknown argument type.")
471 a_inf=util.inf(arg)
472 loc=util.length(arg-a_inf).internal_minGlobalDataPoint() # This gives us the location but not coords
473 x=arg.getFunctionSpace().getX()
474 x_min=x.getTupleForGlobalDataPoint(*loc)
475 return Locator(arg.getFunctionSpace(),x_min)
476
477 def getSupLocator(arg):
478 """
479 Return a Locator for a point with the sup value over all arg.
480 """
481 if not isinstance(arg, escore.Data):
482 raise TypeError("getSupLocator: Unknown argument type.")
483 a_inf=util.sup(arg)
484 loc=util.length(arg-a_inf).internal_minGlobalDataPoint() # This gives us the location but not coords
485 x=arg.getFunctionSpace().getX()
486 x_min=x.getTupleForGlobalDataPoint(*loc)
487 return Locator(arg.getFunctionSpace(),x_min)
488
489
490 class SolverSchemeException(Exception):
491 """
492 This is a generic exception thrown by solvers.
493 """
494 pass
495
496 class IndefinitePreconditioner(SolverSchemeException):
497 """
498 Exception thrown if the preconditioner is not positive definite.
499 """
500 pass
501
502 class MaxIterReached(SolverSchemeException):
503 """
504 Exception thrown if the maximum number of iteration steps is reached.
505 """
506 pass
507
508 class CorrectionFailed(SolverSchemeException):
509 """
510 Exception thrown if no convergence has been achieved in the solution
511 correction scheme.
512 """
513 pass
514
515 class IterationBreakDown(SolverSchemeException):
516 """
517 Exception thrown if the iteration scheme encountered an incurable breakdown.
518 """
519 pass
520
521 class NegativeNorm(SolverSchemeException):
522 """
523 Exception thrown if a norm calculation returns a negative norm.
524 """
525 pass
526
527 def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
528 """
529 Solver for
530
531 *Ax=b*
532
533 with a symmetric and positive definite operator A (more details required!).
534 It uses the conjugate gradient method with preconditioner M providing an
535 approximation of A.
536
537 The iteration is terminated if
538
539 *|r| <= atol+rtol*|r0|*
540
541 where *r0* is the initial residual and *|.|* is the energy norm. In fact
542
543 *|r| = sqrt( bilinearform(Msolve(r),r))*
544
545 For details on the preconditioned conjugate gradient method see the book:
546
547 "Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
548 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
549 C. Romine, and H. van der Vorst".
550
551 :param r: initial residual *r=b-Ax*. ``r`` is altered.
552 :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
553 :param x: an initial guess for the solution
554 :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
555 :param Aprod: returns the value Ax
556 :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
557 argument ``x``. The returned object needs to be of the same type
558 like argument ``r``.
559 :param Msolve: solves Mx=r
560 :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
561 argument ``r``. The returned object needs to be of the same
562 type like argument ``x``.
563 :param bilinearform: inner product ``<x,r>``
564 :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
565 type like argument ``x`` and ``r`` is. The returned value
566 is a ``float``.
567 :param atol: absolute tolerance
568 :type atol: non-negative ``float``
569 :param rtol: relative tolerance
570 :type rtol: non-negative ``float``
571 :param iter_max: maximum number of iteration steps
572 :type iter_max: ``int``
573 :return: the solution approximation and the corresponding residual
574 :rtype: ``tuple``
575 :warning: ``r`` and ``x`` are altered.
576 """
577 iter=0
578 rhat=Msolve(r)
579 d = rhat
580 rhat_dot_r = bilinearform(rhat, r)
581 if rhat_dot_r<0: raise NegativeNorm("negative norm.")
582 norm_r0=math.sqrt(rhat_dot_r)
583 atol2=atol+rtol*norm_r0
584 if atol2<=0:
585 raise ValueError("Non-positive tolarance.")
586 atol2=max(atol2, 100. * util.EPSILON * norm_r0)
587
588 if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
589
590
591 while not math.sqrt(rhat_dot_r) <= atol2:
592 iter+=1
593 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
594
595 q=Aprod(d)
596 alpha = rhat_dot_r / bilinearform(d, q)
597 x += alpha * d
598 if isinstance(q,ArithmeticTuple):
599 r += q * (-alpha) # Doing it the other way calls the float64.__mul__ not AT.__rmul__
600 else:
601 r += (-alpha) * q
602 rhat=Msolve(r)
603 rhat_dot_r_new = bilinearform(rhat, r)
604 beta = rhat_dot_r_new / rhat_dot_r
605 rhat+=beta * d
606 d=rhat
607
608 rhat_dot_r = rhat_dot_r_new
609 if rhat_dot_r<0: raise NegativeNorm("negative norm.")
610 if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
611 if verbose: print(("PCG: tolerance reached after %s steps."%iter))
612 return x,r,math.sqrt(rhat_dot_r)
613
614 class Defect(object):
615 """
616 Defines a non-linear defect F(x) of a variable x. This class includes
617 two functions (bilinearform and eval) that must be overridden by subclassing
618 before use.
619 """
620 def __init__(self):
621 """
622 Initializes defect.
623 """
624 self.setDerivativeIncrementLength()
625
626 def bilinearform(self, x0, x1):
627 """
628 Returns the inner product of x0 and x1
629
630 NOTE: MUST BE OVERRIDDEN BY A SUBCLASS
631
632 :param x0: value for x0
633 :param x1: value for x1
634 :return: the inner product of x0 and x1
635 :rtype: ``float``
636 """
637 raise NotImplementedError("Defect bilinearform method not overridden")
638
639 def norm(self,x):
640 """
641 Returns the norm of argument ``x``.
642
643 :param x: a value
644 :return: norm of argument x
645 :rtype: ``float``
646 :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
647 """
648 s=self.bilinearform(x,x)
649 if s<0: raise NegativeNorm("negative norm.")
650 return math.sqrt(s)
651
652 def eval(self,x):
653 """
654 Returns the value F of a given ``x``.
655
656 NOTE: MUST BE OVERRIDDEN BY A SUBCLASS
657
658 :param x: value for which the defect ``F`` is evaluated
659 :return: value of the defect at ``x``
660 """
661 raise NotImplementedError("Defect eval() method not overridden")
662
663 def __call__(self,x):
664 return self.eval(x)
665
666 def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
667 """
668 Sets the relative length of the increment used to approximate the
669 derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
670 direction of v with x as a starting point.
671
672 :param inc: relative increment length
673 :type inc: positive ``float``
674 """
675 if inc<=0: raise ValueError("positive increment required.")
676 self.__inc=inc
677
678 def getDerivativeIncrementLength(self):
679 """
680 Returns the relative increment length used to approximate the
681 derivative of the defect.
682 :return: value of the defect at ``x``
683 :rtype: positive ``float``
684 """
685 return self.__inc
686
687 def derivative(self, F0, x0, v, v_is_normalised=True):
688 """
689 Returns the directional derivative at ``x0`` in the direction of ``v``.
690
691 :param F0: value of this defect at x0
692 :param x0: value at which derivative is calculated
693 :param v: direction
694 :param v_is_normalised: True to indicate that ``v`` is nomalized
695 (self.norm(v)=0)
696 :return: derivative of this defect at x0 in the direction of ``v``
697 :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
698 used but this method maybe overwritten to use exact evaluation.
699 """
700 normx=self.norm(x0)
701 if normx>0:
702 epsnew = self.getDerivativeIncrementLength() * normx
703 else:
704 epsnew = self.getDerivativeIncrementLength()
705 if not v_is_normalised:
706 normv=self.norm(v)
707 if normv<=0:
708 return F0*0
709 else:
710 epsnew /= normv
711 F1=self.eval(x0 + epsnew * v)
712 return (F1-F0)/epsnew
713
714 ######################################
715 def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
716 """
717 Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
718 criterion:
719
720 *norm(F(x) <= atol + rtol * norm(F(x0)*
721
722 where *x0* is the initial guess.
723
724 :param defect: object defining the function *F*. ``defect.norm`` defines the
725 *norm* used in the stopping criterion.
726 :type defect: `Defect`
727 :param x: initial guess for the solution, ``x`` is altered.
728 :type x: any object type allowing basic operations such as
729 ``numpy.ndarray``, `Data`
730 :param iter_max: maximum number of iteration steps
731 :type iter_max: positive ``int``
732 :param sub_iter_max: maximum number of inner iteration steps
733 :type sub_iter_max: positive ``int``
734 :param atol: absolute tolerance for the solution
735 :type atol: positive ``float``
736 :param rtol: relative tolerance for the solution
737 :type rtol: positive ``float``
738 :param gamma: tolerance safety factor for inner iteration
739 :type gamma: positive ``float``, less than 1
740 :param subtol_max: upper bound for inner tolerance
741 :type subtol_max: positive ``float``, less than 1
742 :return: an approximation of the solution with the desired accuracy
743 :rtype: same type as the initial guess
744 """
745 lmaxit=iter_max
746 if atol<0: raise ValueError("atol needs to be non-negative.")
747 if rtol<0: raise ValueError("rtol needs to be non-negative.")
748 if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
749 if gamma<=0 or gamma>=1: raise ValueError("tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma)
750 if subtol_max<=0 or subtol_max>=1: raise ValueError("upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max)
751
752 F=defect(x)
753 fnrm=defect.norm(F)
754 stop_tol=atol + rtol*fnrm
755 subtol=subtol_max
756 if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
757 if verbose: print((" tolerance = %e."%subtol))
758 iter=1
759 #
760 # main iteration loop
761 #
762 while not fnrm<=stop_tol:
763 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
764 #
765 # adjust subtol_
766 #
767 if iter > 1:
768 rat=fnrm/fnrmo
769 subtol_old=subtol
770 subtol=gamma*rat**2
771 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
772 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
773 #
774 # calculate newton increment xc
775 # if iter_max in __FDGMRES is reached MaxIterReached is thrown
776 # if iter_restart -1 is returned as sub_iter
777 # if atol is reached sub_iter returns the numer of steps performed to get there
778 #
779 #
780 if verbose: print((" subiteration (GMRES) is called with relative tolerance %e."%subtol))
781 try:
782 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
783 except MaxIterReached:
784 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
785 if sub_iter<0:
786 iter+=sub_iter_max
787 else:
788 iter+=sub_iter
789 # ====
790 x+=xc
791 F=defect(x)
792 iter+=1
793 fnrmo, fnrm=fnrm, defect.norm(F)
794 if verbose: print((" step %s: residual %e."%(iter,fnrm)))
795 if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
796 return x
797
798 def __givapp(c,s,vin):
799 """
800 Applies a sequence of Givens rotations (c,s) recursively to the vector
801 ``vin``
802
803 :warning: ``vin`` is altered.
804 """
805 vrot=vin
806 if isinstance(c,float):
807 vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
808 else:
809 for i in range(len(c)):
810 w1=c[i]*vrot[i]-s[i]*vrot[i+1]
811 w2=s[i]*vrot[i]+c[i]*vrot[i+1]
812 vrot[i]=w1
813 vrot[i+1]=w2
814 return vrot
815
816 def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
817 h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
818 c=numpy.zeros(iter_restart,numpy.float64)
819 s=numpy.zeros(iter_restart,numpy.float64)
820 g=numpy.zeros(iter_restart,numpy.float64)
821 v=[]
822
823 rho=defect.norm(F0)
824 if rho<=0.: return x0*0
825
826 v.append(-F0/rho)
827 g[0]=rho
828 iter=0
829 while rho > atol and iter<iter_restart-1:
830 if iter >= iter_max:
831 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
832
833 p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
834 v.append(p)
835
836 v_norm1=defect.norm(v[iter+1])
837
838 # Modified Gram-Schmidt
839 for j in range(iter+1):
840 h[j,iter]=defect.bilinearform(v[j],v[iter+1])
841 v[iter+1]-=h[j,iter]*v[j]
842
843 h[iter+1,iter]=defect.norm(v[iter+1])
844 v_norm2=h[iter+1,iter]
845
846 # Reorthogonalize if needed
847 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
848 for j in range(iter+1):
849 hr=defect.bilinearform(v[j],v[iter+1])
850 h[j,iter]=h[j,iter]+hr
851 v[iter+1] -= hr*v[j]
852
853 v_norm2=defect.norm(v[iter+1])
854 h[iter+1,iter]=v_norm2
855 # watch out for happy breakdown
856 if not v_norm2 == 0:
857 v[iter+1]=v[iter+1]/h[iter+1,iter]
858
859 # Form and store the information for the new Givens rotation
860 if iter > 0 :
861 hhat=numpy.zeros(iter+1,numpy.float64)
862 for i in range(iter+1) : hhat[i]=h[i,iter]
863 hhat=__givapp(c[0:iter],s[0:iter],hhat);
864 for i in range(iter+1) : h[i,iter]=hhat[i]
865
866 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
867
868 if mu!=0 :
869 c[iter]=h[iter,iter]/mu
870 s[iter]=-h[iter+1,iter]/mu
871 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
872 h[iter+1,iter]=0.0
873 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
874 g[iter]=gg[0]
875 g[iter+1]=gg[1]
876
877 # Update the residual norm
878 rho=abs(g[iter+1])
879 iter+=1
880
881 # At this point either iter > iter_max or rho < tol.
882 # It's time to compute x and leave.
883 if iter > 0 :
884 y=numpy.zeros(iter,numpy.float64)
885 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
886 if iter > 1 :
887 i=iter-2
888 while i>=0 :
889 y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
890 i=i-1
891 xhat=v[iter-1]*y[iter-1]
892 for i in range(iter-1):
893 xhat += v[i]*y[i]
894 else :
895 xhat=v[0] * 0
896
897 if iter<iter_restart-1:
898 stopped=iter
899 else:
900 stopped=-1
901
902 return xhat,stopped
903
904 def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
905 """
906 Solver for
907
908 *Ax=b*
909
910 with a general operator A (more details required!).
911 It uses the generalized minimum residual method (GMRES).
912
913 The iteration is terminated if
914
915 *|r| <= atol+rtol*|r0|*
916
917 where *r0* is the initial residual and *|.|* is the energy norm. In fact
918
919 *|r| = sqrt( bilinearform(r,r))*
920
921 :param r: initial residual *r=b-Ax*. ``r`` is altered.
922 :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
923 :param x: an initial guess for the solution
924 :type x: same like ``r``
925 :param Aprod: returns the value Ax
926 :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
927 argument ``x``. The returned object needs to be of the same
928 type like argument ``r``.
929 :param bilinearform: inner product ``<x,r>``
930 :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
931 type like argument ``x`` and ``r``. The returned value is
932 a ``float``.
933 :param atol: absolute tolerance
934 :type atol: non-negative ``float``
935 :param rtol: relative tolerance
936 :type rtol: non-negative ``float``
937 :param iter_max: maximum number of iteration steps
938 :type iter_max: ``int``
939 :param iter_restart: in order to save memory the orthogonalization process
940 is terminated after ``iter_restart`` steps and the
941 iteration is restarted.
942 :type iter_restart: ``int``
943 :return: the solution approximation and the corresponding residual
944 :rtype: ``tuple``
945 :warning: ``r`` and ``x`` are altered.
946 """
947 m=iter_restart
948 restarted=False
949 iter=0
950 if rtol>0:
951 r_dot_r = bilinearform(r, r)
952 if r_dot_r<0: raise NegativeNorm("negative norm.")
953 atol2=atol+rtol*math.sqrt(r_dot_r)
954 if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
955 else:
956 atol2=atol
957 if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
958 if atol2<=0:
959 raise ValueError("Non-positive tolarance.")
960
961 while True:
962 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
963 if restarted:
964 r2 = r-Aprod(x-x2)
965 else:
966 r2=1*r
967 x2=x*1.
968 x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
969 iter+=iter_restart
970 if stopped: break
971 if verbose: print("GMRES: restart.")
972 restarted=True
973 if verbose: print("GMRES: tolerance has been reached.")
974 return x
975
976 def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
977 iter=0
978
979 h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
980 c=numpy.zeros(iter_restart,numpy.float64)
981 s=numpy.zeros(iter_restart,numpy.float64)
982 g=numpy.zeros(iter_restart+1,numpy.float64)
983 v=[]
984
985 r_dot_r = bilinearform(r, r)
986 if r_dot_r<0: raise NegativeNorm("negative norm.")
987 rho=math.sqrt(r_dot_r)
988
989 v.append(r/rho)
990 g[0]=rho
991
992 if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
993 while not (rho<=atol or iter==iter_restart):
994
995 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
996
997 if P_R!=None:
998 p=Aprod(P_R(v[iter]))
999 else:
1000 p=Aprod(v[iter])
1001 v.append(p)
1002
1003 v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1004
1005 # Modified Gram-Schmidt
1006 for j in range(iter+1):
1007 h[j,iter]=bilinearform(v[j],v[iter+1])
1008 v[iter+1]-=h[j,iter]*v[j]
1009
1010 h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
1011 v_norm2=h[iter+1,iter]
1012
1013 # Reorthogonalize if needed
1014 if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default)
1015 for j in range(iter+1):
1016 hr=bilinearform(v[j],v[iter+1])
1017 h[j,iter]=h[j,iter]+hr
1018 v[iter+1] -= hr*v[j]
1019
1020 v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1021 h[iter+1,iter]=v_norm2
1022
1023 # watch out for happy breakdown
1024 if not v_norm2 == 0:
1025 v[iter+1]=v[iter+1]/h[iter+1,iter]
1026
1027 # Form and store the information for the new Givens rotation
1028 if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1029 mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1030
1031 if mu!=0 :
1032 c[iter]=h[iter,iter]/mu
1033 s[iter]=-h[iter+1,iter]/mu
1034 h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1035 h[iter+1,iter]=0.0
1036 gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1037 g[iter]=gg[0]
1038 g[iter+1]=gg[1]
1039 # Update the residual norm
1040
1041 rho=abs(g[iter+1])
1042 if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1043 iter+=1
1044
1045 # At this point either iter > iter_max or rho < tol.
1046 # It's time to compute x and leave.
1047
1048 if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1049 if iter > 0 :
1050 y=numpy.zeros(iter,numpy.float64)
1051 y[iter-1] = g[iter-1] / h[iter-1,iter-1]
1052 if iter > 1 :
1053 i=iter-2
1054 while i>=0 :
1055 y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
1056 i=i-1
1057 xhat=v[iter-1]*y[iter-1]
1058 for i in range(iter-1):
1059 xhat += v[i]*y[i]
1060 else:
1061 xhat=v[0] * 0
1062 if P_R!=None:
1063 x += P_R(xhat)
1064 else:
1065 x += xhat
1066 if iter<iter_restart-1:
1067 stopped=True
1068 else:
1069 stopped=False
1070
1071 return x,stopped
1072
1073 def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1074 """
1075 Solver for
1076
1077 *Ax=b*
1078
1079 with a symmetric and positive definite operator A (more details required!).
1080 It uses the minimum residual method (MINRES) with preconditioner M
1081 providing an approximation of A.
1082
1083 The iteration is terminated if
1084
1085 *|r| <= atol+rtol*|r0|*
1086
1087 where *r0* is the initial residual and *|.|* is the energy norm. In fact
1088
1089 *|r| = sqrt( bilinearform(Msolve(r),r))*
1090
1091 For details on the preconditioned conjugate gradient method see the book:
1092
1093 "Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1094 T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1095 C. Romine, and H. van der Vorst".
1096
1097 :param r: initial residual *r=b-Ax*. ``r`` is altered.
1098 :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1099 :param x: an initial guess for the solution
1100 :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1101 :param Aprod: returns the value Ax
1102 :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1103 argument ``x``. The returned object needs to be of the same
1104 type like argument ``r``.
1105 :param Msolve: solves Mx=r
1106 :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1107 argument ``r``. The returned object needs to be of the same
1108 type like argument ``x``.
1109 :param bilinearform: inner product ``<x,r>``
1110 :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1111 type like argument ``x`` and ``r`` is. The returned value
1112 is a ``float``.
1113 :param atol: absolute tolerance
1114 :type atol: non-negative ``float``
1115 :param rtol: relative tolerance
1116 :type rtol: non-negative ``float``
1117 :param iter_max: maximum number of iteration steps
1118 :type iter_max: ``int``
1119 :return: the solution approximation and the corresponding residual
1120 :rtype: ``tuple``
1121 :warning: ``r`` and ``x`` are altered.
1122 """
1123 #------------------------------------------------------------------
1124 # Set up y and v for the first Lanczos vector v1.
1125 # y = beta1 P' v1, where P = C**(-1).
1126 # v is really P' v1.
1127 #------------------------------------------------------------------
1128 r1 = r
1129 y = Msolve(r)
1130 beta1 = bilinearform(y,r)
1131
1132 if beta1< 0: raise NegativeNorm("negative norm.")
1133
1134 # If r = 0 exactly, stop with x
1135 if beta1==0: return x
1136
1137 if beta1> 0: beta1 = math.sqrt(beta1)
1138
1139 #------------------------------------------------------------------
1140 # Initialize quantities.
1141 # ------------------------------------------------------------------
1142 iter = 0
1143 Anorm = 0
1144 ynorm = 0
1145 oldb = 0
1146 beta = beta1
1147 dbar = 0
1148 epsln = 0
1149 phibar = beta1
1150 rhs1 = beta1
1151 rhs2 = 0
1152 rnorm = phibar
1153 tnorm2 = 0
1154 ynorm2 = 0
1155 cs = -1
1156 sn = 0
1157 w = r*0.
1158 w2 = r*0.
1159 r2 = r1
1160 eps = 0.0001
1161
1162 #---------------------------------------------------------------------
1163 # Main iteration loop.
1164 # --------------------------------------------------------------------
1165 while not rnorm<=atol+rtol*Anorm*ynorm: # checks ||r|| < (||A|| ||x||) * TOL
1166
1167 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1168 iter = iter + 1
1169
1170 #-----------------------------------------------------------------
1171 # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1172 # The general iteration is similar to the case k = 1 with v0 = 0:
1173 #
1174 # p1 = Operator * v1 - beta1 * v0,
1175 # alpha1 = v1'p1,
1176 # q2 = p2 - alpha1 * v1,
1177 # beta2^2 = q2'q2,
1178 # v2 = (1/beta2) q2.
1179 #
1180 # Again, y = betak P vk, where P = C**(-1).
1181 #-----------------------------------------------------------------
1182 s = 1/beta # Normalize previous vector (in y).
1183 v = s*y # v = vk if P = I
1184
1185 y = Aprod(v)
1186
1187 if iter >= 2:
1188 y = y - (beta/oldb)*r1
1189
1190 alfa = bilinearform(v,y) # alphak
1191 y += (- alfa/beta)*r2
1192 r1 = r2
1193 r2 = y
1194 y = Msolve(r2)
1195 oldb = beta # oldb = betak
1196 beta = bilinearform(y,r2) # beta = betak+1^2
1197 if beta < 0: raise NegativeNorm("negative norm.")
1198
1199 beta = math.sqrt( beta )
1200 tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1201
1202 if iter==1: # Initialize a few things.
1203 gmax = abs( alfa ) # alpha1
1204 gmin = gmax # alpha1
1205
1206 # Apply previous rotation Qk-1 to get
1207 # [deltak epslnk+1] = [cs sn][dbark 0 ]
1208 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
1209
1210 oldeps = epsln
1211 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
1212 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
1213 epsln = sn * beta # epsln2 = 0 epslnk+1
1214 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
1215
1216 # Compute the next plane rotation Qk
1217
1218 gamma = math.sqrt(gbar*gbar+beta*beta) # gammak
1219 gamma = max(gamma,eps)
1220 cs = gbar / gamma # ck
1221 sn = beta / gamma # sk
1222 phi = cs * phibar # phik
1223 phibar = sn * phibar # phibark+1
1224
1225 # Update x.
1226
1227 denom = 1/gamma
1228 w1 = w2
1229 w2 = w
1230 w = (v - oldeps*w1 - delta*w2) * denom
1231 x += phi*w
1232
1233 # Go round again.
1234
1235 gmax = max(gmax,gamma)
1236 gmin = min(gmin,gamma)
1237 z = rhs1 / gamma
1238 ynorm2 = z*z + ynorm2
1239 rhs1 = rhs2 - delta*z
1240 rhs2 = - epsln*z
1241
1242 # Estimate various norms and test for convergence.
1243
1244 Anorm = math.sqrt( tnorm2 )
1245 ynorm = math.sqrt( ynorm2 )
1246
1247 rnorm = phibar
1248
1249 return x
1250
1251 def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1252 """
1253 Solver for
1254
1255 *Ax=b*
1256
1257 with a general operator A (more details required!).
1258 It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1259
1260 The iteration is terminated if
1261
1262 *|r| <= atol+rtol*|r0|*
1263
1264 where *r0* is the initial residual and *|.|* is the energy norm. In fact
1265
1266 *|r| = sqrt( bilinearform(r,r))*
1267
1268 :param r: initial residual *r=b-Ax*. ``r`` is altered.
1269 :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1270 :param x: an initial guess for the solution
1271 :type x: same like ``r``
1272 :param Aprod: returns the value Ax
1273 :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1274 argument ``x``. The returned object needs to be of the same type
1275 like argument ``r``.
1276 :param bilinearform: inner product ``<x,r>``
1277 :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1278 type like argument ``x`` and ``r``. The returned value is
1279 a ``float``.
1280 :param atol: absolute tolerance
1281 :type atol: non-negative ``float``
1282 :param rtol: relative tolerance
1283 :type rtol: non-negative ``float``
1284 :param iter_max: maximum number of iteration steps
1285 :type iter_max: ``int``
1286 :rtype: ``tuple``
1287 :warning: ``r`` and ``x`` are altered.
1288 """
1289 u1=0
1290 u2=0
1291 y1=0
1292 y2=0
1293
1294 w = r
1295 y1 = r
1296 iter = 0
1297 d = 0
1298 v = Aprod(y1)
1299 u1 = v
1300
1301 theta = 0.0;
1302 eta = 0.0;
1303 rho=bilinearform(r,r)
1304 if rho < 0: raise NegativeNorm("negative norm.")
1305 tau = math.sqrt(rho)
1306 norm_r0=tau
1307 while tau>atol+rtol*norm_r0:
1308 if iter >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1309
1310 sigma = bilinearform(r,v)
1311 if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1312
1313 alpha = rho / sigma
1314
1315 for j in range(2):
1316 #
1317 # Compute y2 and u2 only if you have to
1318 #
1319 if ( j == 1 ):
1320 y2 = y1 - alpha * v
1321 u2 = Aprod(y2)
1322
1323 m = 2 * (iter+1) - 2 + (j+1)
1324 if j==0:
1325 w = w - alpha * u1
1326 d = y1 + ( theta * theta * eta / alpha ) * d
1327 if j==1:
1328 w = w - alpha * u2
1329 d = y2 + ( theta * theta * eta / alpha ) * d
1330
1331 theta = math.sqrt(bilinearform(w,w))/ tau
1332 c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1333 tau = tau * theta * c
1334 eta = c * c * alpha
1335 x = x + eta * d
1336 #
1337 # Try to terminate the iteration at each pass through the loop
1338 #
1339 if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1340
1341 rhon = bilinearform(r,w)
1342 beta = rhon / rho;
1343 rho = rhon;
1344 y1 = w + beta * y2;
1345 u1 = Aprod(y1)
1346 v = u1 + beta * ( u2 + beta * v )
1347
1348 iter += 1
1349
1350 return x
1351
1352
1353 #############################################
1354
1355 class ArithmeticTuple(object):
1356 """
1357 Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1358 ArithmeticTuple and ``a`` is a float.
1359
1360 Example of usage::
1361
1362 from esys.escript import Data
1363 from numpy import array
1364 a=eData(...)
1365 b=array([1.,4.])
1366 x=ArithmeticTuple(a,b)
1367 y=5.*x
1368
1369 """
1370 def __init__(self,*args):
1371 """
1372 Initializes object with elements ``args``.
1373
1374 :param args: tuple of objects that support inplace add (x+=y) and
1375 scaling (x=a*y)
1376 """
1377 self.__items=list(args)
1378
1379 def __len__(self):
1380 """
1381 Returns the number of items.
1382
1383 :return: number of items
1384 :rtype: ``int``
1385 """
1386 return len(self.__items)
1387
1388 def __getitem__(self,index):
1389 """
1390 Returns item at specified position.
1391
1392 :param index: index of item to be returned
1393 :type index: ``int``
1394 :return: item with index ``index``
1395 """
1396 return self.__items.__getitem__(index)
1397
1398 def __mul__(self,other):
1399 """
1400 Scales by ``other`` from the right.
1401
1402 :param other: scaling factor
1403 :type other: ``float``
1404 :return: itemwise self*other
1405 :rtype: `ArithmeticTuple`
1406 """
1407 out=[]
1408 try:
1409 l=len(other)
1410 if l!=len(self):
1411 raise ValueError("length of arguments don't match.")
1412 for i in range(l):
1413 if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1414 out.append(escore.Data())
1415 else:
1416 out.append(self[i]*other[i])
1417 except TypeError:
1418 for i in range(len(self)):
1419 if self.__isEmpty(self[i]) or self.__isEmpty(other):
1420 out.append(escore.Data())
1421 else:
1422 out.append(self[i]*other)
1423 return ArithmeticTuple(*tuple(out))
1424
1425 def __rmul__(self,other):
1426 """
1427 Scales by ``other`` from the left.
1428
1429 :param other: scaling factor
1430 :type other: ``float``
1431 :return: itemwise other*self
1432 :rtype: `ArithmeticTuple`
1433 """
1434 out=[]
1435 try:
1436 l=len(other)
1437 if l!=len(self):
1438 raise ValueError("length of arguments don't match.")
1439 for i in range(l):
1440 if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1441 out.append(escore.Data())
1442 else:
1443 out.append(other[i]*self[i])
1444 except TypeError:
1445 for i in range(len(self)):
1446 if self.__isEmpty(self[i]) or self.__isEmpty(other):
1447 out.append(escore.Data())
1448 else:
1449 out.append(other*self[i])
1450 return ArithmeticTuple(*tuple(out))
1451
1452 def __div__(self,other):
1453 """
1454 Scales by (1/``other``) from the right.
1455
1456 :param other: scaling factor
1457 :type other: ``float``
1458 :return: itemwise self/other
1459 :rtype: `ArithmeticTuple`
1460 """
1461 return self*(1/other)
1462
1463 def __rdiv__(self,other):
1464 """
1465 Scales by (1/``other``) from the left.
1466
1467 :param other: scaling factor
1468 :type other: ``float``
1469 :return: itemwise other/self
1470 :rtype: `ArithmeticTuple`
1471 """
1472 out=[]
1473 try:
1474 l=len(other)
1475 if l!=len(self):
1476 raise ValueError("length of arguments don't match.")
1477
1478 for i in range(l):
1479 if self.__isEmpty(self[i]):
1480 raise ZeroDivisionError("in component %s"%i)
1481 else:
1482 if self.__isEmpty(other[i]):
1483 out.append(escore.Data())
1484 else:
1485 out.append(other[i]/self[i])
1486 except TypeError:
1487 for i in range(len(self)):
1488 if self.__isEmpty(self[i]):
1489 raise ZeroDivisionError("in component %s"%i)
1490 else:
1491 if self.__isEmpty(other):
1492 out.append(escore.Data())
1493 else:
1494 out.append(other/self[i])
1495 return ArithmeticTuple(*tuple(out))
1496
1497 def __iadd__(self,other):
1498 """
1499 Inplace addition of ``other`` to self.
1500
1501 :param other: increment
1502 :type other: ``ArithmeticTuple``
1503 """
1504 if len(self) != len(other):
1505 raise ValueError("tuple lengths must match.")
1506 for i in range(len(self)):
1507 if self.__isEmpty(self.__items[i]):
1508 self.__items[i]=other[i]
1509 else:
1510 self.__items[i]+=other[i]
1511
1512 return self
1513
1514 def __add__(self,other):
1515 """
1516 Adds ``other`` to self.
1517
1518 :param other: increment
1519 :type other: ``ArithmeticTuple``
1520 """
1521 out=[]
1522 try:
1523 l=len(other)
1524 if l!=len(self):
1525 raise ValueError("length of arguments don't match.")
1526 for i in range(l):
1527 if self.__isEmpty(self[i]):
1528 out.append(other[i])
1529 elif self.__isEmpty(other[i]):
1530 out.append(self[i])
1531 else:
1532 out.append(self[i]+other[i])
1533 except TypeError:
1534 for i in range(len(self)):
1535 if self.__isEmpty(self[i]):
1536 out.append(other)
1537 elif self.__isEmpty(other):
1538 out.append(self[i])
1539 else:
1540 out.append(self[i]+other)
1541 return ArithmeticTuple(*tuple(out))
1542
1543 def __sub__(self,other):
1544 """
1545 Subtracts ``other`` from self.
1546
1547 :param other: decrement
1548 :type other: ``ArithmeticTuple``
1549 """
1550 out=[]
1551 try:
1552 l=len(other)
1553 if l!=len(self):
1554 raise ValueError("length of arguments don't match.")
1555 for i in range(l):
1556 if self.__isEmpty(other[i]):
1557 out.append(self[i])
1558 elif self.__isEmpty(self[i]):
1559 out.append(-other[i])
1560 else:
1561 out.append(self[i]-other[i])
1562 except TypeError:
1563 for i in range(len(self)):
1564 if self.__isEmpty(other):
1565 out.append(self[i])
1566 elif self.__isEmpty(self[i]):
1567 out.append(-other)
1568 else:
1569 out.append(self[i]-other)
1570
1571 return ArithmeticTuple(*tuple(out))
1572
1573 def __isub__(self,other):
1574 """
1575 Inplace subtraction of ``other`` from self.
1576
1577 :param other: decrement
1578 :type other: ``ArithmeticTuple``
1579 """
1580 if len(self) != len(other):
1581 raise ValueError("tuple length must match.")
1582 for i in range(len(self)):
1583 if not self.__isEmpty(other[i]):
1584 if self.__isEmpty(self.__items[i]):
1585 self.__items[i]=-other[i]
1586 else:
1587 self.__items[i]=other[i]
1588 return self
1589
1590 def __neg__(self):
1591 """
1592 Negates values.
1593 """
1594 out=[]
1595 for i in range(len(self)):
1596 if self.__isEmpty(self[i]):
1597 out.append(escore.Data())
1598 else:
1599 out.append(-self[i])
1600
1601 return ArithmeticTuple(*tuple(out))
1602 def __isEmpty(self, d):
1603 if isinstance(d, escore.Data):
1604 return d.isEmpty()
1605 else:
1606 return False
1607
1608 def __str__(self):
1609 s="("
1610 for i in self:
1611 s=s+str(i)+", "
1612 s=s+")"
1613 return s
1614
1615 class HomogeneousSaddlePointProblem(object):
1616 """
1617 This class provides a framework for solving linear homogeneous saddle
1618 point problems of the form::
1619
1620 *Av+B^*p=f*
1621 *Bv =0*
1622
1623 for the unknowns *v* and *p* and given operators *A* and *B* and
1624 given right hand side *f*. *B^** is the adjoint operator of *B*.
1625 *A* may depend weakly on *v* and *p*.
1626 """
1627 def __init__(self, **kwargs):
1628 """
1629 initializes the saddle point problem
1630 """
1631 self.resetControlParameters()
1632 self.setTolerance()
1633 self.setAbsoluteTolerance()
1634 def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1635 """
1636 sets a control parameter
1637
1638 :param K_p: initial value for constant to adjust pressure tolerance
1639 :type K_p: ``float``
1640 :param K_v: initial value for constant to adjust velocity tolerance
1641 :type K_v: ``float``
1642 :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1643 :type rtol_max: ``float``
1644 :param chi_max: maximum tolerable converegence rate.
1645 :type chi_max: ``float``
1646 :param reduction_factor: reduction factor for adjustment factors.
1647 :type reduction_factor: ``float``
1648 """
1649 self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1650
1651 def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1652 """
1653 sets a control parameter
1654
1655
1656 :param K_p: initial value for constant to adjust pressure tolerance
1657 :type K_p: ``float``
1658 :param K_v: initial value for constant to adjust velocity tolerance
1659 :type K_v: ``float``
1660 :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1661 :type rtol_max: ``float``
1662 :param chi_max: maximum tolerable converegence rate.
1663 :type chi_max: ``float``
1664 :type reduction_factor: ``float``
1665 """
1666 if not K_p is None:
1667 if K_p<1:
1668 raise ValueError("K_p need to be greater or equal to 1.")
1669 else:
1670 K_p=self.__K_p
1671
1672 if not K_v is None:
1673 if K_v<1:
1674 raise ValueError("K_v need to be greater or equal to 1.")
1675 else:
1676 K_v=self.__K_v
1677
1678 if not rtol_max is None:
1679 if rtol_max<=0 or rtol_max>=1:
1680 raise ValueError("rtol_max needs to be positive and less than 1.")
1681 else:
1682 rtol_max=self.__rtol_max
1683
1684 if not rtol_min is None:
1685 if rtol_min<=0 or rtol_min>=1:
1686 raise ValueError("rtol_min needs to be positive and less than 1.")
1687 else:
1688 rtol_min=self.__rtol_min
1689
1690 if not chi_max is None:
1691 if chi_max<=0 or chi_max>=1:
1692 raise ValueError("chi_max needs to be positive and less than 1.")
1693 else:
1694 chi_max = self.__chi_max
1695
1696 if not reduction_factor is None:
1697 if reduction_factor<=0 or reduction_factor>1:
1698 raise ValueError("reduction_factor need to be between zero and one.")
1699 else:
1700 reduction_factor=self.__reduction_factor
1701
1702 if not theta is None:
1703 if theta<=0 or theta>1:
1704 raise ValueError("theta need to be between zero and one.")
1705 else:
1706 theta=self.__theta
1707
1708 if rtol_min>=rtol_max:
1709 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1710 self.__chi_max = chi_max
1711 self.__rtol_max = rtol_max
1712 self.__K_p = K_p
1713 self.__K_v = K_v
1714 self.__reduction_factor = reduction_factor
1715 self.__theta = theta
1716 self.__rtol_min=rtol_min
1717
1718 #=============================================================
1719 def inner_pBv(self,p,Bv):
1720 """
1721 Returns inner product of element p and Bv (overwrite).
1722
1723 :param p: a pressure increment
1724 :param Bv: a residual
1725 :return: inner product of element p and Bv
1726 :rtype: ``float``
1727 :note: used if PCG is applied.
1728 """
1729 raise NotImplementedError("no inner product for p and Bv implemented.")
1730
1731 def inner_p(self,p0,p1):
1732 """
1733 Returns inner product of p0 and p1 (overwrite).
1734
1735 :param p0: a pressure
1736 :param p1: a pressure
1737 :return: inner product of p0 and p1
1738 :rtype: ``float``
1739 """
1740 raise NotImplementedError("no inner product for p implemented.")
1741
1742 def norm_v(self,v):
1743 """
1744 Returns the norm of v (overwrite).
1745
1746 :param v: a velovity
1747 :return: norm of v
1748 :rtype: non-negative ``float``
1749 """
1750 raise NotImplementedError("no norm of v implemented.")
1751 def getDV(self, p, v, tol):
1752 """
1753 return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1754
1755 :param p: pressure
1756 :param v: pressure
1757 :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1758 :note: Only *A* may depend on *v* and *p*
1759 """
1760 raise NotImplementedError("no dv calculation implemented.")
1761
1762
1763 def Bv(self,v, tol):
1764 """
1765 Returns Bv with accuracy `tol` (overwrite)
1766
1767 :rtype: equal to the type of p
1768 :note: boundary conditions on p should be zero!
1769 """
1770 raise NotImplementedError("no operator B implemented.")
1771
1772 def norm_Bv(self,Bv):
1773 """
1774 Returns the norm of Bv (overwrite).
1775
1776 :rtype: equal to the type of p
1777 :note: boundary conditions on p should be zero!
1778 """
1779 raise NotImplementedError("no norm of Bv implemented.")
1780
1781 def solve_AinvBt(self,dp, tol):
1782 """
1783 Solves *A dv=B^*dp* with accuracy `tol`
1784
1785 :param dp: a pressure increment
1786 :return: the solution of *A dv=B^*dp*
1787 :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1788 """
1789 raise NotImplementedError("no operator A implemented.")
1790
1791 def solve_prec(self,Bv, tol):
1792 """
1793 Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1794
1795 :rtype: equal to the type of p
1796 :note: boundary conditions on p should be zero!
1797 """
1798 raise NotImplementedError("no preconditioner for Schur complement implemented.")
1799 #=============================================================
1800 def __Aprod_PCG(self,dp):
1801 dv=self.solve_AinvBt(dp, self.__subtol)
1802 return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1803
1804 def __inner_PCG(self,p,r):
1805 return self.inner_pBv(p,r[1])
1806
1807 def __Msolve_PCG(self,r):
1808 return self.solve_prec(r[1], self.__subtol)
1809 #=============================================================
1810 def __Aprod_GMRES(self,p):
1811 return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1812 def __inner_GMRES(self,p0,p1):
1813 return self.inner_p(p0,p1)
1814
1815 #=============================================================
1816 def norm_p(self,p):
1817 """
1818 calculates the norm of ``p``
1819
1820 :param p: a pressure
1821 :return: the norm of ``p`` using the inner product for pressure
1822 :rtype: ``float``
1823 """
1824 f=self.inner_p(p,p)
1825 if f<0: raise ValueError("negative pressure norm.")
1826 return math.sqrt(f)
1827
1828 def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1829 """
1830 Solves the saddle point problem using initial guesses v and p.
1831
1832 :param v: initial guess for velocity
1833 :param p: initial guess for pressure
1834 :type v: `Data`
1835 :type p: `Data`
1836 :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1837 :param max_iter: maximum number of iteration steps per correction
1838 attempt
1839 :param verbose: if True, shows information on the progress of the
1840 saddlepoint problem solver.
1841 :param iter_restart: restart the iteration after ``iter_restart`` steps
1842 (only used if useUzaw=False)
1843 :type usePCG: ``bool``
1844 :type max_iter: ``int``
1845 :type verbose: ``bool``
1846 :type iter_restart: ``int``
1847 :rtype: ``tuple`` of `Data` objects
1848 :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1849 """
1850 return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1851
1852 def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1853 """
1854 see `_solve` method.
1855 """
1856 self.verbose=verbose
1857 rtol=self.getTolerance()
1858 atol=self.getAbsoluteTolerance()
1859
1860 K_p=self.__K_p
1861 K_v=self.__K_v
1862 correction_step=0
1863 converged=False
1864 chi=None
1865 eps=None
1866
1867 if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1868 while not converged:
1869
1870 # get tolerance for velecity increment:
1871 if chi is None:
1872 rtol_v=self.__rtol_max
1873 else:
1874 rtol_v=min(chi/K_v,self.__rtol_max)
1875 rtol_v=max(rtol_v, self.__rtol_min)
1876 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1877 # get velocity increment:
1878 dv1=self.getDV(p,v,rtol_v)
1879 v1=v+dv1
1880 Bv1=self.Bv(v1, rtol_v)
1881 norm_Bv1=self.norm_Bv(Bv1)
1882 norm_dv1=self.norm_v(dv1)
1883 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1884 if norm_dv1*self.__theta < norm_Bv1:
1885 # get tolerance for pressure increment:
1886 large_Bv1=True
1887 if chi is None or eps is None:
1888 rtol_p=self.__rtol_max
1889 else:
1890 rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1891 self.__subtol=max(rtol_p**2, self.__rtol_min)
1892 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1893 # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1894 if usePCG:
1895 dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1896 v2=r[0]
1897 Bv2=r[1]
1898 else:
1899 # don't use!!!!
1900 dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1901 dv2=self.solve_AinvBt(dp, self.__subtol)
1902 v2=v1-dv2
1903 Bv2=self.Bv(v2, self.__subtol)
1904 p2=p+dp
1905 else:
1906 large_Bv1=False
1907 v2=v1
1908 p2=p
1909 # update business:
1910 norm_dv2=self.norm_v(v2-v)
1911 norm_v2=self.norm_v(v2)
1912 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1913 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1914 if eps_old is None:
1915 chi, chi_old = None, chi
1916 else:
1917 chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1918 if eps != None:
1919 if chi !=None:
1920 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1921 else:
1922 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1923 if eps <= rtol*norm_v2+atol :
1924 converged = True
1925 else:
1926 if correction_step>=max_correction_steps:
1927 raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1928 if chi_old!=None:
1929 K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1930 K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1931 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1932 correction_step+=1
1933 v,p =v2, p2
1934 if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1935 return v,p
1936 #========================================================================
1937 def setTolerance(self,tolerance=1.e-4):
1938 """
1939 Sets the relative tolerance for (v,p).
1940
1941 :param tolerance: tolerance to be used
1942 :type tolerance: non-negative ``float``
1943 """
1944 if tolerance<0:
1945 raise ValueError("tolerance must be positive.")
1946 self.__rtol=tolerance
1947
1948 def getTolerance(self):
1949 """
1950 Returns the relative tolerance.
1951
1952 :return: relative tolerance
1953 :rtype: ``float``
1954 """
1955 return self.__rtol
1956
1957 def setAbsoluteTolerance(self,tolerance=0.):
1958 """
1959 Sets the absolute tolerance.
1960
1961 :param tolerance: tolerance to be used
1962 :type tolerance: non-negative ``float``
1963 """
1964 if tolerance<0:
1965 raise ValueError("tolerance must be non-negative.")
1966 self.__atol=tolerance
1967
1968 def getAbsoluteTolerance(self):
1969 """
1970 Returns the absolute tolerance.
1971
1972 :return: absolute tolerance
1973 :rtype: ``float``
1974 """
1975 return self.__atol
1976
1977
1978 def MaskFromBoundaryTag(domain,*tags):
1979 """
1980 Creates a mask on the Solution(domain) function space where the value is
1981 one for samples that touch the boundary tagged by tags.
1982
1983 Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1984
1985 :param domain: domain to be used
1986 :type domain: `escript.Domain`
1987 :param tags: boundary tags
1988 :type tags: ``str``
1989 :return: a mask which marks samples that are touching the boundary tagged
1990 by any of the given tags
1991 :rtype: `escript.Data` of rank 0
1992 """
1993 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1994 d=escore.Scalar(0.,escore.FunctionOnBoundary(domain))
1995 for t in tags: d.setTaggedValue(t,1.)
1996 pde.setValue(y=d)
1997 return util.whereNonZero(pde.getRightHandSide())
1998
1999 def MaskFromTag(domain,*tags):
2000 """
2001 Creates a mask on the Solution(domain) function space where the value is
2002 one for samples that touch regions tagged by tags.
2003
2004 Usage: m=MaskFromTag(domain, "ham")
2005
2006 :param domain: domain to be used
2007 :type domain: `escript.Domain`
2008 :param tags: boundary tags
2009 :type tags: ``str``
2010 :return: a mask which marks samples that are touching the boundary tagged
2011 by any of the given tags
2012 :rtype: `escript.Data` of rank 0
2013 """
2014 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
2015 d=escore.Scalar(0.,escore.Function(domain))
2016 for t in tags: d.setTaggedValue(t,1.)
2017 pde.setValue(Y=d)
2018 return util.whereNonZero(pde.getRightHandSide())
2019
2020 def BoundaryValuesFromVolumeTag(domain,**values):
2021 """
2022 Creates a mask on the Solution(domain) function space where the value is
2023 one for samples that touch regions tagged by tags.
2024
2025 Usage: m=BoundaryValuesFromVolumeTag(domain, ham=1, f=6)
2026
2027 :param domain: domain to be used
2028 :type domain: `escript.Domain`
2029 :return: a mask which marks samples that are touching the boundary tagged
2030 by any of the given tags
2031 :rtype: `escript.Data` of rank 0
2032 """
2033 pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
2034 out=escore.Scalar(0.,escore.FunctionOnBoundary(domain))
2035 for t,v in values.items():
2036 d=escore.Scalar(0.,escore.Function(domain))
2037 d.setTaggedValue(t,1.)
2038 pde.setValue(Y=d)
2039 out+=v*util.whereZero(util.interpolate(util.whereNonZero(pde.getRightHandSide()), escore.FunctionOnBoundary(domain))-1.)
2040 return out
2041

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26