1 |
# |
2 |
# $Id$ |
3 |
# |
4 |
####################################################### |
5 |
# |
6 |
# Copyright 2003-2007 by ACceSS MNRF |
7 |
# Copyright 2007 by University of Queensland |
8 |
# |
9 |
# http://esscc.uq.edu.au |
10 |
# Primary Business: Queensland, Australia |
11 |
# Licensed under the Open Software License version 3.0 |
12 |
# http://www.opensource.org/licenses/osl-3.0.php |
13 |
# |
14 |
####################################################### |
15 |
# |
16 |
|
17 |
""" |
18 |
Provides some tools related to PDEs. |
19 |
|
20 |
Currently includes: |
21 |
- Projector - to project a discontinuous |
22 |
- Locator - to trace values in data objects at a certain location |
23 |
- TimeIntegrationManager - to handel extraplotion in time |
24 |
- SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme |
25 |
|
26 |
@var __author__: name of author |
27 |
@var __copyright__: copyrights |
28 |
@var __license__: licence agreement |
29 |
@var __url__: url entry point on documentation |
30 |
@var __version__: version |
31 |
@var __date__: date of the version |
32 |
""" |
33 |
|
34 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
35 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
36 |
http://www.access.edu.au |
37 |
Primary Business: Queensland, Australia""" |
38 |
__license__="""Licensed under the Open Software License version 3.0 |
39 |
http://www.opensource.org/licenses/osl-3.0.php""" |
40 |
__url__="http://www.iservo.edu.au/esys" |
41 |
__version__="$Revision$" |
42 |
__date__="$Date$" |
43 |
|
44 |
|
45 |
import escript |
46 |
import linearPDEs |
47 |
import numarray |
48 |
import util |
49 |
import math |
50 |
|
51 |
##### Added by Artak |
52 |
# from Numeric import zeros,Int,Float64 |
53 |
################################### |
54 |
|
55 |
|
56 |
class TimeIntegrationManager: |
57 |
""" |
58 |
a simple mechanism to manage time dependend values. |
59 |
|
60 |
typical usage is:: |
61 |
|
62 |
dt=0.1 # time increment |
63 |
tm=TimeIntegrationManager(inital_value,p=1) |
64 |
while t<1. |
65 |
v_guess=tm.extrapolate(dt) # extrapolate to t+dt |
66 |
v=... |
67 |
tm.checkin(dt,v) |
68 |
t+=dt |
69 |
|
70 |
@note: currently only p=1 is supported. |
71 |
""" |
72 |
def __init__(self,*inital_values,**kwargs): |
73 |
""" |
74 |
sets up the value manager where inital_value is the initial value and p is order used for extrapolation |
75 |
""" |
76 |
if kwargs.has_key("p"): |
77 |
self.__p=kwargs["p"] |
78 |
else: |
79 |
self.__p=1 |
80 |
if kwargs.has_key("time"): |
81 |
self.__t=kwargs["time"] |
82 |
else: |
83 |
self.__t=0. |
84 |
self.__v_mem=[inital_values] |
85 |
self.__order=0 |
86 |
self.__dt_mem=[] |
87 |
self.__num_val=len(inital_values) |
88 |
|
89 |
def getTime(self): |
90 |
return self.__t |
91 |
def getValue(self): |
92 |
out=self.__v_mem[0] |
93 |
if len(out)==1: |
94 |
return out[0] |
95 |
else: |
96 |
return out |
97 |
|
98 |
def checkin(self,dt,*values): |
99 |
""" |
100 |
adds new values to the manager. the p+1 last value get lost |
101 |
""" |
102 |
o=min(self.__order+1,self.__p) |
103 |
self.__order=min(self.__order+1,self.__p) |
104 |
v_mem_new=[values] |
105 |
dt_mem_new=[dt] |
106 |
for i in range(o-1): |
107 |
v_mem_new.append(self.__v_mem[i]) |
108 |
dt_mem_new.append(self.__dt_mem[i]) |
109 |
v_mem_new.append(self.__v_mem[o-1]) |
110 |
self.__order=o |
111 |
self.__v_mem=v_mem_new |
112 |
self.__dt_mem=dt_mem_new |
113 |
self.__t+=dt |
114 |
|
115 |
def extrapolate(self,dt): |
116 |
""" |
117 |
extrapolates to dt forward in time. |
118 |
""" |
119 |
if self.__order==0: |
120 |
out=self.__v_mem[0] |
121 |
else: |
122 |
out=[] |
123 |
for i in range(self.__num_val): |
124 |
out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i]) |
125 |
|
126 |
if len(out)==0: |
127 |
return None |
128 |
elif len(out)==1: |
129 |
return out[0] |
130 |
else: |
131 |
return out |
132 |
|
133 |
|
134 |
class Projector: |
135 |
""" |
136 |
The Projector is a factory which projects a discontiuous function onto a |
137 |
continuous function on the a given domain. |
138 |
""" |
139 |
def __init__(self, domain, reduce = True, fast=True): |
140 |
""" |
141 |
Create a continuous function space projector for a domain. |
142 |
|
143 |
@param domain: Domain of the projection. |
144 |
@param reduce: Flag to reduce projection order (default is True) |
145 |
@param fast: Flag to use a fast method based on matrix lumping (default is true) |
146 |
""" |
147 |
self.__pde = linearPDEs.LinearPDE(domain) |
148 |
if fast: |
149 |
self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING) |
150 |
self.__pde.setSymmetryOn() |
151 |
self.__pde.setReducedOrderTo(reduce) |
152 |
self.__pde.setValue(D = 1.) |
153 |
return |
154 |
|
155 |
def __call__(self, input_data): |
156 |
""" |
157 |
Projects input_data onto a continuous function |
158 |
|
159 |
@param input_data: The input_data to be projected. |
160 |
""" |
161 |
out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution()) |
162 |
self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data()) |
163 |
if input_data.getRank()==0: |
164 |
self.__pde.setValue(Y = input_data) |
165 |
out=self.__pde.getSolution() |
166 |
elif input_data.getRank()==1: |
167 |
for i0 in range(input_data.getShape()[0]): |
168 |
self.__pde.setValue(Y = input_data[i0]) |
169 |
out[i0]=self.__pde.getSolution() |
170 |
elif input_data.getRank()==2: |
171 |
for i0 in range(input_data.getShape()[0]): |
172 |
for i1 in range(input_data.getShape()[1]): |
173 |
self.__pde.setValue(Y = input_data[i0,i1]) |
174 |
out[i0,i1]=self.__pde.getSolution() |
175 |
elif input_data.getRank()==3: |
176 |
for i0 in range(input_data.getShape()[0]): |
177 |
for i1 in range(input_data.getShape()[1]): |
178 |
for i2 in range(input_data.getShape()[2]): |
179 |
self.__pde.setValue(Y = input_data[i0,i1,i2]) |
180 |
out[i0,i1,i2]=self.__pde.getSolution() |
181 |
else: |
182 |
for i0 in range(input_data.getShape()[0]): |
183 |
for i1 in range(input_data.getShape()[1]): |
184 |
for i2 in range(input_data.getShape()[2]): |
185 |
for i3 in range(input_data.getShape()[3]): |
186 |
self.__pde.setValue(Y = input_data[i0,i1,i2,i3]) |
187 |
out[i0,i1,i2,i3]=self.__pde.getSolution() |
188 |
return out |
189 |
|
190 |
class NoPDE: |
191 |
""" |
192 |
solves the following problem for u: |
193 |
|
194 |
M{kronecker[i,j]*D[j]*u[j]=Y[i]} |
195 |
|
196 |
with constraint |
197 |
|
198 |
M{u[j]=r[j]} where M{q[j]>0} |
199 |
|
200 |
where D, Y, r and q are given functions of rank 1. |
201 |
|
202 |
In the case of scalars this takes the form |
203 |
|
204 |
M{D*u=Y} |
205 |
|
206 |
with constraint |
207 |
|
208 |
M{u=r} where M{q>0} |
209 |
|
210 |
where D, Y, r and q are given scalar functions. |
211 |
|
212 |
The constraint is overwriting any other condition. |
213 |
|
214 |
@note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention |
215 |
that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole |
216 |
thing is a bit strange and I blame Robert.Woodcock@csiro.au for this. |
217 |
""" |
218 |
def __init__(self,domain,D=None,Y=None,q=None,r=None): |
219 |
""" |
220 |
initialize the problem |
221 |
|
222 |
@param domain: domain of the PDE. |
223 |
@type domain: L{Domain} |
224 |
@param D: coefficient of the solution. |
225 |
@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
226 |
@param Y: right hand side |
227 |
@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
228 |
@param q: location of constraints |
229 |
@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
230 |
@param r: value of solution at locations of constraints |
231 |
@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
232 |
""" |
233 |
self.__domain=domain |
234 |
self.__D=D |
235 |
self.__Y=Y |
236 |
self.__q=q |
237 |
self.__r=r |
238 |
self.__u=None |
239 |
self.__function_space=escript.Solution(self.__domain) |
240 |
def setReducedOn(self): |
241 |
""" |
242 |
sets the L{FunctionSpace} of the solution to L{ReducedSolution} |
243 |
""" |
244 |
self.__function_space=escript.ReducedSolution(self.__domain) |
245 |
self.__u=None |
246 |
|
247 |
def setReducedOff(self): |
248 |
""" |
249 |
sets the L{FunctionSpace} of the solution to L{Solution} |
250 |
""" |
251 |
self.__function_space=escript.Solution(self.__domain) |
252 |
self.__u=None |
253 |
|
254 |
def setValue(self,D=None,Y=None,q=None,r=None): |
255 |
""" |
256 |
assigns values to the parameters. |
257 |
|
258 |
@param D: coefficient of the solution. |
259 |
@type D: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
260 |
@param Y: right hand side |
261 |
@type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
262 |
@param q: location of constraints |
263 |
@type q: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
264 |
@param r: value of solution at locations of constraints |
265 |
@type r: C{float}, C{int}, L{numarray.NumArray}, L{Data} |
266 |
""" |
267 |
if not D==None: |
268 |
self.__D=D |
269 |
self.__u=None |
270 |
if not Y==None: |
271 |
self.__Y=Y |
272 |
self.__u=None |
273 |
if not q==None: |
274 |
self.__q=q |
275 |
self.__u=None |
276 |
if not r==None: |
277 |
self.__r=r |
278 |
self.__u=None |
279 |
|
280 |
def getSolution(self): |
281 |
""" |
282 |
returns the solution |
283 |
|
284 |
@return: the solution of the problem |
285 |
@rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}. |
286 |
""" |
287 |
if self.__u==None: |
288 |
if self.__D==None: |
289 |
raise ValueError,"coefficient D is undefined" |
290 |
D=escript.Data(self.__D,self.__function_space) |
291 |
if D.getRank()>1: |
292 |
raise ValueError,"coefficient D must have rank 0 or 1" |
293 |
if self.__Y==None: |
294 |
self.__u=escript.Data(0.,D.getShape(),self.__function_space) |
295 |
else: |
296 |
self.__u=util.quotient(self.__Y,D) |
297 |
if not self.__q==None: |
298 |
q=util.wherePositive(escript.Data(self.__q,self.__function_space)) |
299 |
self.__u*=(1.-q) |
300 |
if not self.__r==None: self.__u+=q*self.__r |
301 |
return self.__u |
302 |
|
303 |
class Locator: |
304 |
""" |
305 |
Locator provides access to the values of data objects at a given |
306 |
spatial coordinate x. |
307 |
|
308 |
In fact, a Locator object finds the sample in the set of samples of a |
309 |
given function space or domain where which is closest to the given |
310 |
point x. |
311 |
""" |
312 |
|
313 |
def __init__(self,where,x=numarray.zeros((3,))): |
314 |
""" |
315 |
Initializes a Locator to access values in Data objects on the Doamin |
316 |
or FunctionSpace where for the sample point which |
317 |
closest to the given point x. |
318 |
|
319 |
@param where: function space |
320 |
@type where: L{escript.FunctionSpace} |
321 |
@param x: coefficient of the solution. |
322 |
@type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray} |
323 |
""" |
324 |
if isinstance(where,escript.FunctionSpace): |
325 |
self.__function_space=where |
326 |
else: |
327 |
self.__function_space=escript.ContinuousFunction(where) |
328 |
if isinstance(x, list): |
329 |
self.__id=[] |
330 |
for p in x: |
331 |
self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint()) |
332 |
else: |
333 |
self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint() |
334 |
|
335 |
def __str__(self): |
336 |
""" |
337 |
Returns the coordinates of the Locator as a string. |
338 |
""" |
339 |
x=self.getX() |
340 |
if instance(x,list): |
341 |
out="[" |
342 |
first=True |
343 |
for xx in x: |
344 |
if not first: |
345 |
out+="," |
346 |
else: |
347 |
first=False |
348 |
out+=str(xx) |
349 |
out+="]>" |
350 |
else: |
351 |
out=str(x) |
352 |
return out |
353 |
|
354 |
def getX(self): |
355 |
""" |
356 |
Returns the exact coordinates of the Locator. |
357 |
""" |
358 |
return self(self.getFunctionSpace().getX()) |
359 |
|
360 |
def getFunctionSpace(self): |
361 |
""" |
362 |
Returns the function space of the Locator. |
363 |
""" |
364 |
return self.__function_space |
365 |
|
366 |
def getId(self,item=None): |
367 |
""" |
368 |
Returns the identifier of the location. |
369 |
""" |
370 |
if item == None: |
371 |
return self.__id |
372 |
else: |
373 |
if isinstance(self.__id,list): |
374 |
return self.__id[item] |
375 |
else: |
376 |
return self.__id |
377 |
|
378 |
|
379 |
def __call__(self,data): |
380 |
""" |
381 |
Returns the value of data at the Locator of a Data object otherwise |
382 |
the object is returned. |
383 |
""" |
384 |
return self.getValue(data) |
385 |
|
386 |
def getValue(self,data): |
387 |
""" |
388 |
Returns the value of data at the Locator if data is a Data object |
389 |
otherwise the object is returned. |
390 |
""" |
391 |
if isinstance(data,escript.Data): |
392 |
if data.getFunctionSpace()==self.getFunctionSpace(): |
393 |
dat=data |
394 |
else: |
395 |
dat=data.interpolate(self.getFunctionSpace()) |
396 |
id=self.getId() |
397 |
r=data.getRank() |
398 |
if isinstance(id,list): |
399 |
out=[] |
400 |
for i in id: |
401 |
o=data.getValueOfGlobalDataPoint(*i) |
402 |
if data.getRank()==0: |
403 |
out.append(o[0]) |
404 |
else: |
405 |
out.append(o) |
406 |
return out |
407 |
else: |
408 |
out=data.getValueOfGlobalDataPoint(*id) |
409 |
if data.getRank()==0: |
410 |
return out[0] |
411 |
else: |
412 |
return out |
413 |
else: |
414 |
return data |
415 |
|
416 |
class SolverSchemeException(Exception): |
417 |
""" |
418 |
exceptions thrown by solvers |
419 |
""" |
420 |
pass |
421 |
|
422 |
class IndefinitePreconditioner(SolverSchemeException): |
423 |
""" |
424 |
the preconditioner is not positive definite. |
425 |
""" |
426 |
pass |
427 |
class MaxIterReached(SolverSchemeException): |
428 |
""" |
429 |
maxium number of iteration steps is reached. |
430 |
""" |
431 |
pass |
432 |
class IterationBreakDown(SolverSchemeException): |
433 |
""" |
434 |
iteration scheme econouters an incurable breakdown. |
435 |
""" |
436 |
pass |
437 |
class NegativeNorm(SolverSchemeException): |
438 |
""" |
439 |
a norm calculation returns a negative norm. |
440 |
""" |
441 |
pass |
442 |
|
443 |
class IterationHistory(object): |
444 |
""" |
445 |
The IterationHistory class is used to define a stopping criterium. It keeps track of the |
446 |
residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by |
447 |
a given tolerance. |
448 |
""" |
449 |
def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False): |
450 |
""" |
451 |
Initialization |
452 |
|
453 |
@param tolerance: tolerance |
454 |
@type tolerance: positive C{float} |
455 |
@param verbose: switches on the printing out some information |
456 |
@type verbose: C{bool} |
457 |
""" |
458 |
if not tolerance>0.: |
459 |
raise ValueError,"tolerance needs to be positive." |
460 |
self.tolerance=tolerance |
461 |
self.verbose=verbose |
462 |
self.history=[] |
463 |
def stoppingcriterium(self,norm_r,r,x): |
464 |
""" |
465 |
returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]} is the residual norm at the first call. |
466 |
|
467 |
|
468 |
@param norm_r: current residual norm |
469 |
@type norm_r: non-negative C{float} |
470 |
@param r: current residual (not used) |
471 |
@param x: current solution approximation (not used) |
472 |
@return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned. |
473 |
@rtype: C{bool} |
474 |
|
475 |
""" |
476 |
self.history.append(norm_r) |
477 |
if self.verbose: print "iter: %s: inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1]) |
478 |
return self.history[-1]<=self.tolerance * self.history[0] |
479 |
|
480 |
def stoppingcriterium2(self,norm_r,norm_b): |
481 |
""" |
482 |
returns True if the C{norm_r} is C{tolerance}*C{norm_b} |
483 |
|
484 |
|
485 |
@param norm_r: current residual norm |
486 |
@type norm_r: non-negative C{float} |
487 |
@param norm_b: norm of right hand side |
488 |
@type norm_b: non-negative C{float} |
489 |
@return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned. |
490 |
@rtype: C{bool} |
491 |
|
492 |
""" |
493 |
self.history.append(norm_r) |
494 |
if self.verbose: print "iter: %s: norm(r) = %e"%(len(self.history)-1, self.history[-1]) |
495 |
return self.history[-1]<=self.tolerance * norm_b |
496 |
|
497 |
def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100): |
498 |
""" |
499 |
Solver for |
500 |
|
501 |
M{Ax=b} |
502 |
|
503 |
with a symmetric and positive definite operator A (more details required!). |
504 |
It uses the conjugate gradient method with preconditioner M providing an approximation of A. |
505 |
|
506 |
The iteration is terminated if the C{stoppingcriterium} function return C{True}. |
507 |
|
508 |
For details on the preconditioned conjugate gradient method see the book: |
509 |
|
510 |
Templates for the Solution of Linear Systems by R. Barrett, M. Berry, |
511 |
T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, |
512 |
C. Romine, and H. van der Vorst. |
513 |
|
514 |
@param b: the right hand side of the liner system. C{b} is altered. |
515 |
@type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y) |
516 |
@param Aprod: returns the value Ax |
517 |
@type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}. |
518 |
@param Msolve: solves Mx=r |
519 |
@type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same |
520 |
type like argument C{x}. |
521 |
@param bilinearform: inner product C{<x,r>} |
522 |
@type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}. |
523 |
@param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step. |
524 |
@type stoppingcriterium: function that returns C{True} or C{False} |
525 |
@param x: an initial guess for the solution. If no C{x} is given 0*b is used. |
526 |
@type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y) |
527 |
@param iter_max: maximum number of iteration steps. |
528 |
@type iter_max: C{int} |
529 |
@return: the solution approximation and the corresponding residual |
530 |
@rtype: C{tuple} |
531 |
@warning: C{b} and C{x} are altered. |
532 |
""" |
533 |
iter=0 |
534 |
if x==None: |
535 |
x=0*b |
536 |
else: |
537 |
b += (-1)*Aprod(x) |
538 |
r=b |
539 |
rhat=Msolve(r) |
540 |
d = rhat |
541 |
rhat_dot_r = bilinearform(rhat, r) |
542 |
if rhat_dot_r<0: raise NegativeNorm,"negative norm." |
543 |
|
544 |
while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x): |
545 |
iter+=1 |
546 |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
547 |
|
548 |
q=Aprod(d) |
549 |
alpha = rhat_dot_r / bilinearform(d, q) |
550 |
x += alpha * d |
551 |
r += (-alpha) * q |
552 |
|
553 |
rhat=Msolve(r) |
554 |
rhat_dot_r_new = bilinearform(rhat, r) |
555 |
beta = rhat_dot_r_new / rhat_dot_r |
556 |
rhat+=beta * d |
557 |
d=rhat |
558 |
|
559 |
rhat_dot_r = rhat_dot_r_new |
560 |
if rhat_dot_r<0: raise NegativeNorm,"negative norm." |
561 |
|
562 |
return x,r |
563 |
|
564 |
|
565 |
############################ |
566 |
# Added by Artak |
567 |
#################################3 |
568 |
|
569 |
#Apply a sequence of k Givens rotations, used within gmres codes |
570 |
# vrot=givapp(c, s, vin, k) |
571 |
def givapp(c,s,vin): |
572 |
vrot=vin # warning: vin is altered!!!! |
573 |
if isinstance(c,float): |
574 |
vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]] |
575 |
else: |
576 |
for i in range(len(c)): |
577 |
w1=c[i]*vrot[i]-s[i]*vrot[i+1] |
578 |
w2=s[i]*vrot[i]+c[i]*vrot[i+1] |
579 |
vrot[i:i+2]=w1,w2 |
580 |
return vrot |
581 |
|
582 |
def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10): |
583 |
m=iter_restart |
584 |
iter=0 |
585 |
while True: |
586 |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
587 |
x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m) |
588 |
iter+=iter_restart |
589 |
if stopped: break |
590 |
return x |
591 |
|
592 |
def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10): |
593 |
iter=0 |
594 |
r=Msolve(b) |
595 |
r_dot_r = bilinearform(r, r) |
596 |
if r_dot_r<0: raise NegativeNorm,"negative norm." |
597 |
norm_b=math.sqrt(r_dot_r) |
598 |
|
599 |
if x==None: |
600 |
x=0*b |
601 |
else: |
602 |
r=Msolve(b-Aprod(x)) |
603 |
r_dot_r = bilinearform(r, r) |
604 |
if r_dot_r<0: raise NegativeNorm,"negative norm." |
605 |
|
606 |
h=numarray.zeros((iter_restart,iter_restart),numarray.Float64) |
607 |
c=numarray.zeros(iter_restart,numarray.Float64) |
608 |
s=numarray.zeros(iter_restart,numarray.Float64) |
609 |
g=numarray.zeros(iter_restart,numarray.Float64) |
610 |
v=[] |
611 |
|
612 |
rho=math.sqrt(r_dot_r) |
613 |
v.append(r/rho) |
614 |
g[0]=rho |
615 |
|
616 |
while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1): |
617 |
|
618 |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
619 |
|
620 |
|
621 |
p=Msolve(Aprod(v[iter])) |
622 |
|
623 |
v.append(p) |
624 |
|
625 |
v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1])) |
626 |
|
627 |
# Modified Gram-Schmidt |
628 |
for j in range(iter+1): |
629 |
h[j][iter]=bilinearform(v[j],v[iter+1]) |
630 |
v[iter+1]+=(-1.)*h[j][iter]*v[j] |
631 |
|
632 |
h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1])) |
633 |
v_norm2=h[iter+1][iter] |
634 |
|
635 |
|
636 |
# Reorthogonalize if needed |
637 |
if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default) |
638 |
for j in range(iter+1): |
639 |
hr=bilinearform(v[j],v[iter+1]) |
640 |
h[j][iter]=h[j][iter]+hr #vhat |
641 |
v[iter+1] +=(-1.)*hr*v[j] |
642 |
|
643 |
v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1])) |
644 |
h[iter+1][iter]=v_norm2 |
645 |
|
646 |
# watch out for happy breakdown |
647 |
if v_norm2 != 0: |
648 |
v[iter+1]=v[iter+1]/h[iter+1][iter] |
649 |
|
650 |
# Form and store the information for the new Givens rotation |
651 |
if iter > 0 : |
652 |
hhat=[] |
653 |
for i in range(iter+1) : hhat.append(h[i][iter]) |
654 |
hhat=givapp(c[0:iter],s[0:iter],hhat); |
655 |
for i in range(iter+1) : h[i][iter]=hhat[i] |
656 |
|
657 |
mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter]) |
658 |
if mu!=0 : |
659 |
c[iter]=h[iter][iter]/mu |
660 |
s[iter]=-h[iter+1][iter]/mu |
661 |
h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter] |
662 |
h[iter+1][iter]=0.0 |
663 |
g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2]) |
664 |
|
665 |
# Update the residual norm |
666 |
rho=abs(g[iter+1]) |
667 |
iter+=1 |
668 |
|
669 |
# At this point either iter > iter_max or rho < tol. |
670 |
# It's time to compute x and leave. |
671 |
|
672 |
if iter > 0 : |
673 |
y=numarray.zeros(iter,numarray.Float64) |
674 |
y[iter-1] = g[iter-1] / h[iter-1][iter-1] |
675 |
if iter > 1 : |
676 |
i=iter-2 |
677 |
while i>=0 : |
678 |
y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i] |
679 |
i=i-1 |
680 |
xhat=v[iter-1]*y[iter-1] |
681 |
for i in range(iter-1): |
682 |
xhat += v[i]*y[i] |
683 |
else : xhat=v[0] |
684 |
|
685 |
x += xhat |
686 |
if iter!=iter_restart-1: |
687 |
stopped=True |
688 |
else: |
689 |
stopped=False |
690 |
|
691 |
return x,stopped |
692 |
|
693 |
def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100): |
694 |
|
695 |
# |
696 |
# minres solves the system of linear equations Ax = b |
697 |
# where A is a symmetric matrix (possibly indefinite or singular) |
698 |
# and b is a given vector. |
699 |
# |
700 |
# "A" may be a dense or sparse matrix (preferably sparse!) |
701 |
# or the name of a function such that |
702 |
# y = A(x) |
703 |
# returns the product y = Ax for any given vector x. |
704 |
# |
705 |
# "M" defines a positive-definite preconditioner M = C C'. |
706 |
# "M" may be a dense or sparse matrix (preferably sparse!) |
707 |
# or the name of a function such that |
708 |
# solves the system My = x for any given vector x. |
709 |
# |
710 |
# |
711 |
|
712 |
# Initialize |
713 |
|
714 |
iter = 0 |
715 |
Anorm = 0 |
716 |
ynorm = 0 |
717 |
x=x*0 |
718 |
#------------------------------------------------------------------ |
719 |
# Set up y and v for the first Lanczos vector v1. |
720 |
# y = beta1 P' v1, where P = C**(-1). |
721 |
# v is really P' v1. |
722 |
#------------------------------------------------------------------ |
723 |
r1 = b |
724 |
y = Msolve(b) |
725 |
beta1 = bilinearform(b,y) |
726 |
|
727 |
if beta1< 0: raise NegativeNorm,"negative norm." |
728 |
|
729 |
# If b = 0 exactly, stop with x = 0. |
730 |
if beta1==0: return x*0. |
731 |
|
732 |
if beta1> 0: |
733 |
beta1 = math.sqrt(beta1) # Normalize y to get v1 later. |
734 |
|
735 |
#------------------------------------------------------------------ |
736 |
# Initialize other quantities. |
737 |
# ------------------------------------------------------------------ |
738 |
oldb = 0 |
739 |
beta = beta1 |
740 |
dbar = 0 |
741 |
epsln = 0 |
742 |
phibar = beta1 |
743 |
rhs1 = beta1 |
744 |
rhs2 = 0 |
745 |
rnorm = phibar |
746 |
tnorm2 = 0 |
747 |
ynorm2 = 0 |
748 |
cs = -1 |
749 |
sn = 0 |
750 |
w = b*0. |
751 |
w2 = b*0. |
752 |
r2 = r1 |
753 |
eps = 0.0001 |
754 |
|
755 |
#--------------------------------------------------------------------- |
756 |
# Main iteration loop. |
757 |
# -------------------------------------------------------------------- |
758 |
while not stoppingcriterium(rnorm,Anorm*ynorm): # ||r|| / (||A|| ||x||) |
759 |
|
760 |
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
761 |
iter = iter + 1 |
762 |
|
763 |
#----------------------------------------------------------------- |
764 |
# Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,... |
765 |
# The general iteration is similar to the case k = 1 with v0 = 0: |
766 |
# |
767 |
# p1 = Operator * v1 - beta1 * v0, |
768 |
# alpha1 = v1'p1, |
769 |
# q2 = p2 - alpha1 * v1, |
770 |
# beta2^2 = q2'q2, |
771 |
# v2 = (1/beta2) q2. |
772 |
# |
773 |
# Again, y = betak P vk, where P = C**(-1). |
774 |
#----------------------------------------------------------------- |
775 |
s = 1/beta # Normalize previous vector (in y). |
776 |
v = s*y # v = vk if P = I |
777 |
|
778 |
y = Aprod(v) |
779 |
|
780 |
if iter >= 2: |
781 |
y = y - (beta/oldb)*r1 |
782 |
|
783 |
alfa = bilinearform(v,y) # alphak |
784 |
y = (- alfa/beta)*r2 + y |
785 |
r1 = r2 |
786 |
r2 = y |
787 |
y = Msolve(r2) |
788 |
oldb = beta # oldb = betak |
789 |
beta = bilinearform(r2,y) # beta = betak+1^2 |
790 |
if beta < 0: raise NegativeNorm,"negative norm." |
791 |
|
792 |
beta = math.sqrt( beta ) |
793 |
tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta |
794 |
|
795 |
if iter==1: # Initialize a few things. |
796 |
gmax = abs( alfa ) # alpha1 |
797 |
gmin = gmax # alpha1 |
798 |
|
799 |
# Apply previous rotation Qk-1 to get |
800 |
# [deltak epslnk+1] = [cs sn][dbark 0 ] |
801 |
# [gbar k dbar k+1] [sn -cs][alfak betak+1]. |
802 |
|
803 |
oldeps = epsln |
804 |
delta = cs * dbar + sn * alfa # delta1 = 0 deltak |
805 |
gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k |
806 |
epsln = sn * beta # epsln2 = 0 epslnk+1 |
807 |
dbar = - cs * beta # dbar 2 = beta2 dbar k+1 |
808 |
|
809 |
# Compute the next plane rotation Qk |
810 |
|
811 |
gamma = math.sqrt(gbar*gbar+beta*beta) # gammak |
812 |
gamma = max(gamma,eps) |
813 |
cs = gbar / gamma # ck |
814 |
sn = beta / gamma # sk |
815 |
phi = cs * phibar # phik |
816 |
phibar = sn * phibar # phibark+1 |
817 |
|
818 |
# Update x. |
819 |
|
820 |
denom = 1/gamma |
821 |
w1 = w2 |
822 |
w2 = w |
823 |
w = (v - oldeps*w1 - delta*w2) * denom |
824 |
x = x + phi*w |
825 |
|
826 |
# Go round again. |
827 |
|
828 |
gmax = max(gmax,gamma) |
829 |
gmin = min(gmin,gamma) |
830 |
z = rhs1 / gamma |
831 |
ynorm2 = z*z + ynorm2 |
832 |
rhs1 = rhs2 - delta*z |
833 |
rhs2 = - epsln*z |
834 |
|
835 |
# Estimate various norms and test for convergence. |
836 |
|
837 |
Anorm = math.sqrt( tnorm2 ) |
838 |
ynorm = math.sqrt( ynorm2 ) |
839 |
|
840 |
rnorm = phibar |
841 |
|
842 |
# Return final answer. |
843 |
return x |
844 |
|
845 |
############################################# |
846 |
|
847 |
class ArithmeticTuple(object): |
848 |
""" |
849 |
tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float. |
850 |
|
851 |
example of usage: |
852 |
|
853 |
from esys.escript import Data |
854 |
from numarray import array |
855 |
a=Data(...) |
856 |
b=array([1.,4.]) |
857 |
x=ArithmeticTuple(a,b) |
858 |
y=5.*x |
859 |
|
860 |
""" |
861 |
def __init__(self,*args): |
862 |
""" |
863 |
initialize object with elements args. |
864 |
|
865 |
@param args: tuple of object that support implace add (x+=y) and scaling (x=a*y) |
866 |
""" |
867 |
self.__items=list(args) |
868 |
|
869 |
def __len__(self): |
870 |
""" |
871 |
number of items |
872 |
|
873 |
@return: number of items |
874 |
@rtype: C{int} |
875 |
""" |
876 |
return len(self.__items) |
877 |
|
878 |
def __getitem__(self,index): |
879 |
""" |
880 |
get an item |
881 |
|
882 |
@param index: item to be returned |
883 |
@type index: C{int} |
884 |
@return: item with index C{index} |
885 |
""" |
886 |
return self.__items.__getitem__(index) |
887 |
|
888 |
def __mul__(self,other): |
889 |
""" |
890 |
scaling from the right |
891 |
|
892 |
@param other: scaling factor |
893 |
@type other: C{float} |
894 |
@return: itemwise self*other |
895 |
@rtype: L{ArithmeticTuple} |
896 |
""" |
897 |
out=[] |
898 |
for i in range(len(self)): |
899 |
out.append(self[i]*other) |
900 |
return ArithmeticTuple(*tuple(out)) |
901 |
|
902 |
def __rmul__(self,other): |
903 |
""" |
904 |
scaling from the left |
905 |
|
906 |
@param other: scaling factor |
907 |
@type other: C{float} |
908 |
@return: itemwise other*self |
909 |
@rtype: L{ArithmeticTuple} |
910 |
""" |
911 |
out=[] |
912 |
for i in range(len(self)): |
913 |
out.append(other*self[i]) |
914 |
return ArithmeticTuple(*tuple(out)) |
915 |
|
916 |
######################### |
917 |
# Added by Artak |
918 |
######################### |
919 |
def __div__(self,other): |
920 |
""" |
921 |
dividing from the right |
922 |
|
923 |
@param other: scaling factor |
924 |
@type other: C{float} |
925 |
@return: itemwise self/other |
926 |
@rtype: L{ArithmeticTuple} |
927 |
""" |
928 |
out=[] |
929 |
for i in range(len(self)): |
930 |
out.append(self[i]/other) |
931 |
return ArithmeticTuple(*tuple(out)) |
932 |
|
933 |
def __rdiv__(self,other): |
934 |
""" |
935 |
dividing from the left |
936 |
|
937 |
@param other: scaling factor |
938 |
@type other: C{float} |
939 |
@return: itemwise other/self |
940 |
@rtype: L{ArithmeticTuple} |
941 |
""" |
942 |
out=[] |
943 |
for i in range(len(self)): |
944 |
out.append(other/self[i]) |
945 |
return ArithmeticTuple(*tuple(out)) |
946 |
|
947 |
##########################################33 |
948 |
|
949 |
def __iadd__(self,other): |
950 |
""" |
951 |
in-place add of other to self |
952 |
|
953 |
@param other: increment |
954 |
@type other: C{ArithmeticTuple} |
955 |
""" |
956 |
if len(self) != len(other): |
957 |
raise ValueError,"tuple length must match." |
958 |
for i in range(len(self)): |
959 |
self.__items[i]+=other[i] |
960 |
return self |
961 |
|
962 |
class HomogeneousSaddlePointProblem(object): |
963 |
""" |
964 |
This provides a framwork for solving homogeneous saddle point problem of the form |
965 |
|
966 |
Av+B^*p=f |
967 |
Bv =0 |
968 |
|
969 |
for the unknowns v and p and given operators A and B and given right hand side f. |
970 |
B^* is the adjoint operator of B is the given inner product. |
971 |
|
972 |
""" |
973 |
def __init__(self,**kwargs): |
974 |
self.setTolerance() |
975 |
self.setToleranceReductionFactor() |
976 |
|
977 |
def initialize(self): |
978 |
""" |
979 |
initialize the problem (overwrite) |
980 |
""" |
981 |
pass |
982 |
def B(self,v): |
983 |
""" |
984 |
returns Bv (overwrite) |
985 |
@rtype: equal to the type of p |
986 |
|
987 |
@note: boundary conditions on p should be zero! |
988 |
""" |
989 |
pass |
990 |
|
991 |
def inner(self,p0,p1): |
992 |
""" |
993 |
returns inner product of two element p0 and p1 (overwrite) |
994 |
|
995 |
@type p0: equal to the type of p |
996 |
@type p1: equal to the type of p |
997 |
@rtype: C{float} |
998 |
|
999 |
@rtype: equal to the type of p |
1000 |
""" |
1001 |
pass |
1002 |
|
1003 |
def solve_A(self,u,p): |
1004 |
""" |
1005 |
solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite) |
1006 |
|
1007 |
@rtype: equal to the type of v |
1008 |
@note: boundary conditions on v should be zero! |
1009 |
""" |
1010 |
pass |
1011 |
|
1012 |
def solve_prec(self,p): |
1013 |
""" |
1014 |
provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite) |
1015 |
|
1016 |
@rtype: equal to the type of p |
1017 |
""" |
1018 |
pass |
1019 |
|
1020 |
def stoppingcriterium(self,Bv,v,p): |
1021 |
""" |
1022 |
returns a True if iteration is terminated. (overwrite) |
1023 |
|
1024 |
@rtype: C{bool} |
1025 |
""" |
1026 |
pass |
1027 |
|
1028 |
def __inner(self,p,r): |
1029 |
return self.inner(p,r[1]) |
1030 |
|
1031 |
def __inner_p(self,p1,p2): |
1032 |
return self.inner(p1,p2) |
1033 |
|
1034 |
def __stoppingcriterium(self,norm_r,r,p): |
1035 |
return self.stoppingcriterium(r[1],r[0],p) |
1036 |
|
1037 |
def __stoppingcriterium_GMRES(self,norm_r,norm_b): |
1038 |
return self.stoppingcriterium_GMRES(norm_r,norm_b) |
1039 |
|
1040 |
def __stoppingcriterium_MINRES(self,norm_r,norm_Ax): |
1041 |
return self.stoppingcriterium_MINRES(norm_r,norm_Ax) |
1042 |
|
1043 |
|
1044 |
def setTolerance(self,tolerance=1.e-8): |
1045 |
self.__tol=tolerance |
1046 |
def getTolerance(self): |
1047 |
return self.__tol |
1048 |
def setToleranceReductionFactor(self,reduction=0.01): |
1049 |
self.__reduction=reduction |
1050 |
def getSubProblemTolerance(self): |
1051 |
return self.__reduction*self.getTolerance() |
1052 |
|
1053 |
def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'): |
1054 |
""" |
1055 |
solves the saddle point problem using initial guesses v and p. |
1056 |
|
1057 |
@param max_iter: maximum number of iteration steps. |
1058 |
""" |
1059 |
self.verbose=verbose |
1060 |
self.show_details=show_details and self.verbose |
1061 |
|
1062 |
# assume p is known: then v=A^-1(f-B^*p) |
1063 |
# which leads to BA^-1B^*p = BA^-1f |
1064 |
|
1065 |
# Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask) |
1066 |
|
1067 |
|
1068 |
self.__z=v+self.solve_A(v,p*0) |
1069 |
|
1070 |
Bz=self.B(self.__z) |
1071 |
# |
1072 |
# solve BA^-1B^*p = Bz |
1073 |
# |
1074 |
# note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv |
1075 |
# |
1076 |
# with Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask) |
1077 |
# A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask) |
1078 |
# |
1079 |
self.iter=0 |
1080 |
if solver=='GMRES': |
1081 |
if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter |
1082 |
p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.) |
1083 |
# solve Au=f-B^*p |
1084 |
# A(u-v)=f-B^*p-Av |
1085 |
# u=v+(u-v) |
1086 |
u=v+self.solve_A(v,p) |
1087 |
|
1088 |
if solver=='MINRES': |
1089 |
if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter |
1090 |
p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.) |
1091 |
# solve Au=f-B^*p |
1092 |
# A(u-v)=f-B^*p-Av |
1093 |
# u=v+(u-v) |
1094 |
u=v+self.solve_A(v,p) |
1095 |
|
1096 |
if solver=='PCG': |
1097 |
if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter |
1098 |
p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p) |
1099 |
u=r[0] |
1100 |
|
1101 |
print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u) |
1102 |
|
1103 |
return u,p |
1104 |
|
1105 |
def __Msolve(self,r): |
1106 |
return self.solve_prec(r[1]) |
1107 |
|
1108 |
def __Msolve_GMRES(self,r): |
1109 |
return self.solve_prec(r) |
1110 |
|
1111 |
|
1112 |
def __Aprod(self,p): |
1113 |
# return BA^-1B*p |
1114 |
#solve Av =-B^*p as Av =f-Az-B^*p |
1115 |
v=self.solve_A(self.__z,-p) |
1116 |
return ArithmeticTuple(v, self.B(v)) |
1117 |
|
1118 |
def __Aprod_GMRES(self,p): |
1119 |
# return BA^-1B*p |
1120 |
#solve Av =-B^*p as Av =f-Az-B^*p |
1121 |
v=self.solve_A(self.__z,-p) |
1122 |
return self.B(v) |
1123 |
|
1124 |
class SaddlePointProblem(object): |
1125 |
""" |
1126 |
This implements a solver for a saddlepoint problem |
1127 |
|
1128 |
M{f(u,p)=0} |
1129 |
M{g(u)=0} |
1130 |
|
1131 |
for u and p. The problem is solved with an inexact Uszawa scheme for p: |
1132 |
|
1133 |
M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})} |
1134 |
M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})} |
1135 |
|
1136 |
where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of |
1137 |
A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper' |
1138 |
Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays |
1139 |
in fact the role of a preconditioner. |
1140 |
""" |
1141 |
def __init__(self,verbose=False,*args): |
1142 |
""" |
1143 |
initializes the problem |
1144 |
|
1145 |
@param verbose: switches on the printing out some information |
1146 |
@type verbose: C{bool} |
1147 |
@note: this method may be overwritten by a particular saddle point problem |
1148 |
""" |
1149 |
if not isinstance(verbose,bool): |
1150 |
raise TypeError("verbose needs to be of type bool.") |
1151 |
self.__verbose=verbose |
1152 |
self.relaxation=1. |
1153 |
|
1154 |
def trace(self,text): |
1155 |
""" |
1156 |
prints text if verbose has been set |
1157 |
|
1158 |
@param text: a text message |
1159 |
@type text: C{str} |
1160 |
""" |
1161 |
if self.__verbose: print "%s: %s"%(str(self),text) |
1162 |
|
1163 |
def solve_f(self,u,p,tol=1.e-8): |
1164 |
""" |
1165 |
solves |
1166 |
|
1167 |
A_f du = f(u,p) |
1168 |
|
1169 |
with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u. |
1170 |
|
1171 |
@param u: current approximation of u |
1172 |
@type u: L{escript.Data} |
1173 |
@param p: current approximation of p |
1174 |
@type p: L{escript.Data} |
1175 |
@param tol: tolerance expected for du |
1176 |
@type tol: C{float} |
1177 |
@return: increment du |
1178 |
@rtype: L{escript.Data} |
1179 |
@note: this method has to be overwritten by a particular saddle point problem |
1180 |
""" |
1181 |
pass |
1182 |
|
1183 |
def solve_g(self,u,tol=1.e-8): |
1184 |
""" |
1185 |
solves |
1186 |
|
1187 |
Q_g dp = g(u) |
1188 |
|
1189 |
with Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. |
1190 |
|
1191 |
@param u: current approximation of u |
1192 |
@type u: L{escript.Data} |
1193 |
@param tol: tolerance expected for dp |
1194 |
@type tol: C{float} |
1195 |
@return: increment dp |
1196 |
@rtype: L{escript.Data} |
1197 |
@note: this method has to be overwritten by a particular saddle point problem |
1198 |
""" |
1199 |
pass |
1200 |
|
1201 |
def inner(self,p0,p1): |
1202 |
""" |
1203 |
inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1) |
1204 |
@return: inner product of p0 and p1 |
1205 |
@rtype: C{float} |
1206 |
""" |
1207 |
pass |
1208 |
|
1209 |
subiter_max=3 |
1210 |
def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None): |
1211 |
""" |
1212 |
runs the solver |
1213 |
|
1214 |
@param u0: initial guess for C{u} |
1215 |
@type u0: L{esys.escript.Data} |
1216 |
@param p0: initial guess for C{p} |
1217 |
@type p0: L{esys.escript.Data} |
1218 |
@param tolerance: tolerance for relative error in C{u} and C{p} |
1219 |
@type tolerance: positive C{float} |
1220 |
@param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance} |
1221 |
@type tolerance_u: positive C{float} |
1222 |
@param iter_max: maximum number of iteration steps. |
1223 |
@type iter_max: C{int} |
1224 |
@param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the |
1225 |
relaxation factor. If C{accepted_reduction=None} no backtracking is used. |
1226 |
@type accepted_reduction: positive C{float} or C{None} |
1227 |
@param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used. |
1228 |
@type relaxation: C{float} or C{None} |
1229 |
""" |
1230 |
tol=1.e-2 |
1231 |
if tolerance_u==None: tolerance_u=tolerance |
1232 |
if not relaxation==None: self.relaxation=relaxation |
1233 |
if accepted_reduction ==None: |
1234 |
angle_limit=0. |
1235 |
elif accepted_reduction>=1.: |
1236 |
angle_limit=0. |
1237 |
else: |
1238 |
angle_limit=util.sqrt(1-accepted_reduction**2) |
1239 |
self.iter=0 |
1240 |
u=u0 |
1241 |
p=p0 |
1242 |
# |
1243 |
# initialize things: |
1244 |
# |
1245 |
converged=False |
1246 |
# |
1247 |
# start loop: |
1248 |
# |
1249 |
# initial search direction is g |
1250 |
# |
1251 |
while not converged : |
1252 |
if self.iter>iter_max: |
1253 |
raise ArithmeticError("no convergence after %s steps."%self.iter) |
1254 |
f_new=self.solve_f(u,p,tol) |
1255 |
norm_f_new = util.Lsup(f_new) |
1256 |
u_new=u-f_new |
1257 |
g_new=self.solve_g(u_new,tol) |
1258 |
self.iter+=1 |
1259 |
norm_g_new = util.sqrt(self.inner(g_new,g_new)) |
1260 |
if norm_f_new==0. and norm_g_new==0.: return u, p |
1261 |
if self.iter>1 and not accepted_reduction==None: |
1262 |
# |
1263 |
# did we manage to reduce the norm of G? I |
1264 |
# if not we start a backtracking procedure |
1265 |
# |
1266 |
# print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g |
1267 |
if norm_g_new > accepted_reduction * norm_g: |
1268 |
sub_iter=0 |
1269 |
s=self.relaxation |
1270 |
d=g |
1271 |
g_last=g |
1272 |
self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s)) |
1273 |
while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g: |
1274 |
dg= g_new-g_last |
1275 |
norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation) |
1276 |
rad=self.inner(g_new,dg)/self.relaxation |
1277 |
# print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit |
1278 |
# print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g |
1279 |
if abs(rad) < norm_dg*norm_g_new * angle_limit: |
1280 |
if sub_iter>0: self.trace(" no further improvements expected from backtracking.") |
1281 |
break |
1282 |
r=self.relaxation |
1283 |
self.relaxation= - rad/norm_dg**2 |
1284 |
s+=self.relaxation |
1285 |
##### |
1286 |
# a=g_new+self.relaxation*dg/r |
1287 |
# print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation |
1288 |
##### |
1289 |
g_last=g_new |
1290 |
p+=self.relaxation*d |
1291 |
f_new=self.solve_f(u,p,tol) |
1292 |
u_new=u-f_new |
1293 |
g_new=self.solve_g(u_new,tol) |
1294 |
self.iter+=1 |
1295 |
norm_f_new = util.Lsup(f_new) |
1296 |
norm_g_new = util.sqrt(self.inner(g_new,g_new)) |
1297 |
# print " ",sub_iter," new g norm",norm_g_new |
1298 |
self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s)) |
1299 |
# |
1300 |
# can we expect reduction of g? |
1301 |
# |
1302 |
# u_last=u_new |
1303 |
sub_iter+=1 |
1304 |
self.relaxation=s |
1305 |
# |
1306 |
# check for convergence: |
1307 |
# |
1308 |
norm_u_new = util.Lsup(u_new) |
1309 |
p_new=p+self.relaxation*g_new |
1310 |
norm_p_new = util.sqrt(self.inner(p_new,p_new)) |
1311 |
self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation)) |
1312 |
|
1313 |
if self.iter>1: |
1314 |
dg2=g_new-g |
1315 |
df2=f_new-f |
1316 |
norm_dg2=util.sqrt(self.inner(dg2,dg2)) |
1317 |
norm_df2=util.Lsup(df2) |
1318 |
# print norm_g_new, norm_g, norm_dg, norm_p, tolerance |
1319 |
tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new |
1320 |
tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new |
1321 |
if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f: |
1322 |
converged=True |
1323 |
f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new |
1324 |
self.trace("convergence after %s steps."%self.iter) |
1325 |
return u,p |
1326 |
# def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.): |
1327 |
# tol=1.e-2 |
1328 |
# iter=0 |
1329 |
# converged=False |
1330 |
# u=u0*1. |
1331 |
# p=p0*1. |
1332 |
# while not converged and iter<iter_max: |
1333 |
# du=self.solve_f(u,p,tol) |
1334 |
# u-=du |
1335 |
# norm_du=util.Lsup(du) |
1336 |
# norm_u=util.Lsup(u) |
1337 |
# |
1338 |
# dp=self.relaxation*self.solve_g(u,tol) |
1339 |
# p+=dp |
1340 |
# norm_dp=util.sqrt(self.inner(dp,dp)) |
1341 |
# norm_p=util.sqrt(self.inner(p,p)) |
1342 |
# print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p) |
1343 |
# iter+=1 |
1344 |
# |
1345 |
# converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p) |
1346 |
# if converged: |
1347 |
# print "convergence after %s steps."%iter |
1348 |
# else: |
1349 |
# raise ArithmeticError("no convergence after %s steps."%iter) |
1350 |
# |
1351 |
# return u,p |
1352 |
|
1353 |
def MaskFromBoundaryTag(function_space,*tags): |
1354 |
""" |
1355 |
create a mask on the given function space which one for samples |
1356 |
that touch the boundary tagged by tags. |
1357 |
|
1358 |
usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right") |
1359 |
|
1360 |
@param function_space: a given function space |
1361 |
@type function_space: L{escript.FunctionSpace} |
1362 |
@param tags: boundray tags |
1363 |
@type tags: C{str} |
1364 |
@return: a mask which marks samples used by C{function_space} that are touching the |
1365 |
boundary tagged by any of the given tags. |
1366 |
@rtype: L{escript.Data} of rank 0 |
1367 |
""" |
1368 |
pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1) |
1369 |
d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain())) |
1370 |
for t in tags: d.setTaggedValue(t,1.) |
1371 |
pde.setValue(y=d) |
1372 |
out=util.whereNonZero(pde.getRightHandSide()) |
1373 |
if out.getFunctionSpace() == function_space: |
1374 |
return out |
1375 |
else: |
1376 |
return util.whereNonZero(util.interpolate(out,function_space)) |
1377 |
|
1378 |
|
1379 |
|