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

Annotation of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 921 - (hide annotations)
Fri Jan 5 00:54:37 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 21477 byte(s)
I have done some clarification on functions that allow to access individual data point values in a Data object. 
The term "data point number" is always local on a MPI process and referes to the value (data_point_in_sample, sample)
as a single identifyer (data_point_in_sample + sample * number_data_points_per_sample). a "global data point number"
referes to a tuple of a processour id and local data point number.

The function convertToNumArrayFromSampleNo has been removed now and convertToNumArrayFromDPNo renamed to getValueOfDataPoint.
There are two new functions:

   getNumberOfDataPoints
   setValueOfDataPoint

This allows you to do things like:

  in=Data(..)
  out=Data(..)
   for i in xrange(in.getNumberOfDataPoints())
       in_loc=in.getValueOfDataPoint(i)
       out_loc=< some operations on in_loc>
       out.setValueOfDataPoint(i,out_loc)


Also mindp  is renamed to  minGlobalDataPoint and there is a new function getValueOfGlobalDataPoint. While in MPI the functions getNumberOfDataPoints and getValueOfDataPoint are working locally on each process (so the code above is executed in parallel).
the latter allows getting a single value across all processors. 


1 jgs 121 # $Id$
2    
3     """
4 jgs 149 Provides some tools related to PDEs.
5 jgs 121
6 jgs 149 Currently includes:
7     - Projector - to project a discontinuous
8 gross 351 - Locator - to trace values in data objects at a certain location
9     - TimeIntegrationManager - to handel extraplotion in time
10 gross 867 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11 gross 637
12     @var __author__: name of author
13     @var __copyright__: copyrights
14     @var __license__: licence agreement
15     @var __url__: url entry point on documentation
16     @var __version__: version
17     @var __date__: date of the version
18 jgs 121 """
19    
20 gross 637 __author__="Lutz Gross, l.gross@uq.edu.au"
21 elspeth 609 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22     http://www.access.edu.au
23     Primary Business: Queensland, Australia"""
24 elspeth 614 __license__="""Licensed under the Open Software License version 3.0
25     http://www.opensource.org/licenses/osl-3.0.php"""
26 gross 637 __url__="http://www.iservo.edu.au/esys"
27     __version__="$Revision$"
28     __date__="$Date$"
29 elspeth 609
30 gross 637
31 jgs 149 import escript
32     import linearPDEs
33 jgs 121 import numarray
34 jgs 149 import util
35 jgs 121
36 gross 351 class TimeIntegrationManager:
37     """
38     a simple mechanism to manage time dependend values.
39    
40 gross 720 typical usage is::
41 gross 351
42 gross 720 dt=0.1 # time increment
43     tm=TimeIntegrationManager(inital_value,p=1)
44     while t<1.
45     v_guess=tm.extrapolate(dt) # extrapolate to t+dt
46     v=...
47     tm.checkin(dt,v)
48     t+=dt
49 gross 351
50 gross 720 @note: currently only p=1 is supported.
51 gross 351 """
52     def __init__(self,*inital_values,**kwargs):
53     """
54     sets up the value manager where inital_value is the initial value and p is order used for extrapolation
55     """
56     if kwargs.has_key("p"):
57     self.__p=kwargs["p"]
58     else:
59     self.__p=1
60     if kwargs.has_key("time"):
61     self.__t=kwargs["time"]
62     else:
63     self.__t=0.
64     self.__v_mem=[inital_values]
65     self.__order=0
66     self.__dt_mem=[]
67     self.__num_val=len(inital_values)
68    
69     def getTime(self):
70     return self.__t
71 gross 396 def getValue(self):
72 gross 409 out=self.__v_mem[0]
73     if len(out)==1:
74     return out[0]
75     else:
76     return out
77    
78 gross 351 def checkin(self,dt,*values):
79     """
80     adds new values to the manager. the p+1 last value get lost
81     """
82     o=min(self.__order+1,self.__p)
83     self.__order=min(self.__order+1,self.__p)
84     v_mem_new=[values]
85     dt_mem_new=[dt]
86     for i in range(o-1):
87     v_mem_new.append(self.__v_mem[i])
88     dt_mem_new.append(self.__dt_mem[i])
89     v_mem_new.append(self.__v_mem[o-1])
90     self.__order=o
91     self.__v_mem=v_mem_new
92     self.__dt_mem=dt_mem_new
93     self.__t+=dt
94    
95     def extrapolate(self,dt):
96     """
97     extrapolates to dt forward in time.
98     """
99     if self.__order==0:
100     out=self.__v_mem[0]
101     else:
102     out=[]
103     for i in range(self.__num_val):
104     out.append((1.+dt/self.__dt_mem[0])*self.__v_mem[0][i]-dt/self.__dt_mem[0]*self.__v_mem[1][i])
105    
106     if len(out)==0:
107     return None
108     elif len(out)==1:
109     return out[0]
110     else:
111     return out
112 gross 396
113 gross 351
114 jgs 121 class Projector:
115 jgs 149 """
116     The Projector is a factory which projects a discontiuous function onto a
117     continuous function on the a given domain.
118     """
119 jgs 121 def __init__(self, domain, reduce = True, fast=True):
120     """
121 jgs 149 Create a continuous function space projector for a domain.
122 jgs 121
123 jgs 149 @param domain: Domain of the projection.
124     @param reduce: Flag to reduce projection order (default is True)
125     @param fast: Flag to use a fast method based on matrix lumping (default is true)
126 jgs 121 """
127 jgs 149 self.__pde = linearPDEs.LinearPDE(domain)
128 jgs 148 if fast:
129 jgs 149 self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
130 jgs 121 self.__pde.setSymmetryOn()
131     self.__pde.setReducedOrderTo(reduce)
132     self.__pde.setValue(D = 1.)
133     return
134    
135     def __del__(self):
136     return
137    
138     def __call__(self, input_data):
139     """
140 jgs 149 Projects input_data onto a continuous function
141 jgs 121
142 jgs 149 @param input_data: The input_data to be projected.
143 jgs 121 """
144 gross 525 out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
145 jgs 121 if input_data.getRank()==0:
146     self.__pde.setValue(Y = input_data)
147     out=self.__pde.getSolution()
148     elif input_data.getRank()==1:
149     for i0 in range(input_data.getShape()[0]):
150     self.__pde.setValue(Y = input_data[i0])
151     out[i0]=self.__pde.getSolution()
152     elif input_data.getRank()==2:
153     for i0 in range(input_data.getShape()[0]):
154     for i1 in range(input_data.getShape()[1]):
155     self.__pde.setValue(Y = input_data[i0,i1])
156     out[i0,i1]=self.__pde.getSolution()
157     elif input_data.getRank()==3:
158     for i0 in range(input_data.getShape()[0]):
159     for i1 in range(input_data.getShape()[1]):
160     for i2 in range(input_data.getShape()[2]):
161     self.__pde.setValue(Y = input_data[i0,i1,i2])
162     out[i0,i1,i2]=self.__pde.getSolution()
163     else:
164     for i0 in range(input_data.getShape()[0]):
165     for i1 in range(input_data.getShape()[1]):
166     for i2 in range(input_data.getShape()[2]):
167     for i3 in range(input_data.getShape()[3]):
168     self.__pde.setValue(Y = input_data[i0,i1,i2,i3])
169     out[i0,i1,i2,i3]=self.__pde.getSolution()
170     return out
171    
172 gross 525 class NoPDE:
173     """
174     solves the following problem for u:
175 jgs 121
176 gross 525 M{kronecker[i,j]*D[j]*u[j]=Y[i]}
177    
178     with constraint
179    
180     M{u[j]=r[j]} where M{q[j]>0}
181    
182     where D, Y, r and q are given functions of rank 1.
183    
184     In the case of scalars this takes the form
185    
186     M{D*u=Y}
187    
188     with constraint
189    
190     M{u=r} where M{q>0}
191    
192     where D, Y, r and q are given scalar functions.
193    
194     The constraint is overwriting any other condition.
195    
196 gross 720 @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention
197     that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole
198     thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.
199 gross 525 """
200     def __init__(self,domain,D=None,Y=None,q=None,r=None):
201     """
202     initialize the problem
203    
204     @param domain: domain of the PDE.
205     @type domain: L{Domain}
206     @param D: coefficient of the solution.
207 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
208 gross 525 @param Y: right hand side
209 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
210 gross 525 @param q: location of constraints
211 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
212 gross 525 @param r: value of solution at locations of constraints
213 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
214 gross 525 """
215     self.__domain=domain
216     self.__D=D
217     self.__Y=Y
218     self.__q=q
219     self.__r=r
220     self.__u=None
221     self.__function_space=escript.Solution(self.__domain)
222     def setReducedOn(self):
223     """
224     sets the L{FunctionSpace} of the solution to L{ReducedSolution}
225     """
226     self.__function_space=escript.ReducedSolution(self.__domain)
227     self.__u=None
228    
229     def setReducedOff(self):
230     """
231     sets the L{FunctionSpace} of the solution to L{Solution}
232     """
233     self.__function_space=escript.Solution(self.__domain)
234     self.__u=None
235    
236     def setValue(self,D=None,Y=None,q=None,r=None):
237     """
238     assigns values to the parameters.
239    
240     @param D: coefficient of the solution.
241 gross 720 @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}
242 gross 525 @param Y: right hand side
243 gross 720 @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}
244 gross 525 @param q: location of constraints
245 gross 720 @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}
246 gross 525 @param r: value of solution at locations of constraints
247 gross 720 @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}
248 gross 525 """
249     if not D==None:
250     self.__D=D
251     self.__u=None
252     if not Y==None:
253     self.__Y=Y
254     self.__u=None
255     if not q==None:
256     self.__q=q
257     self.__u=None
258     if not r==None:
259     self.__r=r
260     self.__u=None
261    
262     def getSolution(self):
263     """
264     returns the solution
265    
266     @return: the solution of the problem
267     @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.
268     """
269     if self.__u==None:
270     if self.__D==None:
271     raise ValueError,"coefficient D is undefined"
272     D=escript.Data(self.__D,self.__function_space)
273     if D.getRank()>1:
274     raise ValueError,"coefficient D must have rank 0 or 1"
275     if self.__Y==None:
276     self.__u=escript.Data(0.,D.getShape(),self.__function_space)
277     else:
278     self.__u=util.quotient(self.__Y,D)
279     if not self.__q==None:
280     q=util.wherePositive(escript.Data(self.__q,self.__function_space))
281     self.__u*=(1.-q)
282     if not self.__r==None: self.__u+=q*self.__r
283     return self.__u
284    
285 jgs 147 class Locator:
286     """
287 jgs 149 Locator provides access to the values of data objects at a given
288     spatial coordinate x.
289    
290     In fact, a Locator object finds the sample in the set of samples of a
291     given function space or domain where which is closest to the given
292     point x.
293 jgs 147 """
294    
295     def __init__(self,where,x=numarray.zeros((3,))):
296 jgs 149 """
297     Initializes a Locator to access values in Data objects on the Doamin
298     or FunctionSpace where for the sample point which
299     closest to the given point x.
300 gross 880
301     @param where: function space
302     @type where: L{escript.FunctionSpace}
303     @param x: coefficient of the solution.
304     @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}
305 jgs 149 """
306     if isinstance(where,escript.FunctionSpace):
307 jgs 147 self.__function_space=where
308 jgs 121 else:
309 jgs 149 self.__function_space=escript.ContinuousFunction(where)
310 gross 880 if isinstance(x, list):
311     self.__id=[]
312     for p in x:
313 gross 921 self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
314 gross 880 else:
315 gross 921 self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
316 jgs 121
317 jgs 147 def __str__(self):
318 jgs 149 """
319     Returns the coordinates of the Locator as a string.
320     """
321 gross 880 x=self.getX()
322     if instance(x,list):
323     out="["
324     first=True
325     for xx in x:
326     if not first:
327     out+=","
328     else:
329     first=False
330     out+=str(xx)
331     out+="]>"
332     else:
333     out=str(x)
334     return out
335 jgs 121
336 gross 880 def getX(self):
337     """
338     Returns the exact coordinates of the Locator.
339     """
340     return self(self.getFunctionSpace().getX())
341    
342 jgs 147 def getFunctionSpace(self):
343 jgs 149 """
344     Returns the function space of the Locator.
345     """
346 jgs 147 return self.__function_space
347    
348 gross 880 def getId(self,item=None):
349 jgs 149 """
350     Returns the identifier of the location.
351     """
352 gross 880 if item == None:
353     return self.__id
354     else:
355     if isinstance(self.__id,list):
356     return self.__id[item]
357     else:
358     return self.__id
359 jgs 121
360    
361 jgs 147 def __call__(self,data):
362 jgs 149 """
363     Returns the value of data at the Locator of a Data object otherwise
364     the object is returned.
365     """
366 jgs 147 return self.getValue(data)
367 jgs 121
368 jgs 147 def getValue(self,data):
369 jgs 149 """
370     Returns the value of data at the Locator if data is a Data object
371     otherwise the object is returned.
372     """
373     if isinstance(data,escript.Data):
374 jgs 147 if data.getFunctionSpace()==self.getFunctionSpace():
375 gross 880 dat=data
376 jgs 147 else:
377 gross 880 dat=data.interpolate(self.getFunctionSpace())
378     id=self.getId()
379     r=data.getRank()
380     if isinstance(id,list):
381     out=[]
382     for i in id:
383 gross 921 o=data.getValueOfGlobalDataPoint(*i)
384 gross 880 if data.getRank()==0:
385     out.append(o[0])
386     else:
387     out.append(o)
388     return out
389 jgs 147 else:
390 gross 921 out=data.getValueOfGlobalDataPoint(*id)
391 gross 880 if data.getRank()==0:
392     return out[0]
393     else:
394     return out
395 jgs 147 else:
396     return data
397 jgs 149
398 gross 867 class SaddlePointProblem(object):
399     """
400     This implements a solver for a saddlepoint problem
401    
402 gross 877 M{f(u,p)=0}
403     M{g(u)=0}
404 gross 867
405     for u and p. The problem is solved with an inexact Uszawa scheme for p:
406    
407 gross 877 M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})
408     M{Q_g (p^{k+1}-p^{k}) = g(u^{k+1})}
409 gross 867
410     where Q_f is an approximation of the Jacobiean A_f of f with respect to u and Q_f is an approximation of
411     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'
412     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays
413     in fact the role of a preconditioner.
414     """
415     def __init__(self,verbose=False,*args):
416     """
417     initializes the problem
418    
419     @parm verbose: switches on the printing out some information
420     @type verbose: C{bool}
421     @note: this method may be overwritten by a particular saddle point problem
422     """
423     self.__verbose=verbose
424 gross 877 self.relaxation=1.
425 gross 867
426     def trace(self,text):
427     """
428     prints text if verbose has been set
429    
430     @parm text: a text message
431     @type text: C{str}
432     """
433     if self.__verbose: print "%s: %s"%(str(self),text)
434    
435 gross 873 def solve_f(self,u,p,tol=1.e-8):
436 gross 867 """
437     solves
438    
439     A_f du = f(u,p)
440    
441     with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.
442    
443     @param u: current approximation of u
444     @type u: L{escript.Data}
445     @param p: current approximation of p
446     @type p: L{escript.Data}
447 gross 873 @param tol: tolerance expected for du
448 gross 867 @type tol: C{float}
449     @return: increment du
450     @rtype: L{escript.Data}
451     @note: this method has to be overwritten by a particular saddle point problem
452     """
453     pass
454    
455 gross 873 def solve_g(self,u,tol=1.e-8):
456 gross 867 """
457     solves
458    
459     Q_g dp = g(u)
460    
461     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.
462    
463     @param u: current approximation of u
464     @type u: L{escript.Data}
465 gross 873 @param tol: tolerance expected for dp
466     @type tol: C{float}
467 gross 867 @return: increment dp
468     @rtype: L{escript.Data}
469     @note: this method has to be overwritten by a particular saddle point problem
470     """
471     pass
472    
473     def inner(self,p0,p1):
474     """
475     inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)
476     @return: inner product of p0 and p1
477     @rtype: C{float}
478     """
479     pass
480    
481 gross 877 subiter_max=3
482     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
483     """
484     runs the solver
485 gross 873
486 gross 877 @param u0: initial guess for C{u}
487     @type u0: L{esys.escript.Data}
488     @param p0: initial guess for C{p}
489     @type p0: L{esys.escript.Data}
490     @param tolerance: tolerance for relative error in C{u} and C{p}
491     @type tolerance: positive C{float}
492     @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}
493     @type tolerance_u: positive C{float}
494     @param iter_max: maximum number of iteration steps.
495     @type iter_max: C{int}
496     @param accepted_reduction: if the norm g cannot be reduced by C{accepted_reduction} backtracking to adapt the
497     relaxation factor. If C{accepted_reduction=None} no backtracking is used.
498     @type accepted_reduction: positive C{float} or C{None}
499     @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.
500     @type relaxation: C{float} or C{None}
501     """
502     tol=1.e-2
503     if tolerance_u==None: tolerance_u=tolerance
504     if not relaxation==None: self.relaxation=relaxation
505     if accepted_reduction ==None:
506     angle_limit=0.
507     elif accepted_reduction>=1.:
508     angle_limit=0.
509     else:
510     angle_limit=util.sqrt(1-accepted_reduction**2)
511     self.iter=0
512     u=u0
513     p=p0
514     #
515     # initialize things:
516     #
517     converged=False
518     #
519     # start loop:
520     #
521     # initial search direction is g
522     #
523     while not converged :
524     if self.iter>iter_max:
525     raise ArithmeticError("no convergence after %s steps."%self.iter)
526     f_new=self.solve_f(u,p,tol)
527     norm_f_new = util.Lsup(f_new)
528     u_new=u-f_new
529     g_new=self.solve_g(u_new,tol)
530     self.iter+=1
531     norm_g_new = util.sqrt(self.inner(g_new,g_new))
532     if norm_f_new==0. and norm_g_new==0.: return u, p
533     if self.iter>1 and not accepted_reduction==None:
534     #
535     # did we manage to reduce the norm of G? I
536     # if not we start a backtracking procedure
537     #
538     # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
539     if norm_g_new > accepted_reduction * norm_g:
540     sub_iter=0
541     s=self.relaxation
542     d=g
543     g_last=g
544     self.trace(" start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
545     while sub_iter < self.subiter_max and norm_g_new > accepted_reduction * norm_g:
546     dg= g_new-g_last
547     norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
548     rad=self.inner(g_new,dg)/self.relaxation
549     # print " ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
550     # print " ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
551     if abs(rad) < norm_dg*norm_g_new * angle_limit:
552     if sub_iter>0: self.trace(" no further improvements expected from backtracking.")
553     break
554     r=self.relaxation
555     self.relaxation= - rad/norm_dg**2
556     s+=self.relaxation
557     #####
558     # a=g_new+self.relaxation*dg/r
559     # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
560     #####
561     g_last=g_new
562     p+=self.relaxation*d
563     f_new=self.solve_f(u,p,tol)
564     u_new=u-f_new
565     g_new=self.solve_g(u_new,tol)
566     self.iter+=1
567     norm_f_new = util.Lsup(f_new)
568     norm_g_new = util.sqrt(self.inner(g_new,g_new))
569     # print " ",sub_iter," new g norm",norm_g_new
570     self.trace(" %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
571     #
572     # can we expect reduction of g?
573     #
574     # u_last=u_new
575     sub_iter+=1
576     self.relaxation=s
577     #
578     # check for convergence:
579     #
580     norm_u_new = util.Lsup(u_new)
581     p_new=p+self.relaxation*g_new
582     norm_p_new = util.sqrt(self.inner(p_new,p_new))
583     self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))
584 gross 873
585 gross 877 if self.iter>1:
586     dg2=g_new-g
587     df2=f_new-f
588     norm_dg2=util.sqrt(self.inner(dg2,dg2))
589     norm_df2=util.Lsup(df2)
590     # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
591     tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
592     tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
593     if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
594     converged=True
595     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
596     self.trace("convergence after %s steps."%self.iter)
597     return u,p
598     # def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
599     # tol=1.e-2
600     # iter=0
601     # converged=False
602     # u=u0*1.
603     # p=p0*1.
604     # while not converged and iter<iter_max:
605     # du=self.solve_f(u,p,tol)
606     # u-=du
607     # norm_du=util.Lsup(du)
608     # norm_u=util.Lsup(u)
609     #
610     # dp=self.relaxation*self.solve_g(u,tol)
611     # p+=dp
612     # norm_dp=util.sqrt(self.inner(dp,dp))
613     # norm_p=util.sqrt(self.inner(p,p))
614     # 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)
615     # iter+=1
616     #
617     # converged = (norm_du <= tolerance*norm_u) and (norm_dp <= tolerance*norm_p)
618     # if converged:
619     # print "convergence after %s steps."%iter
620     # else:
621     # raise ArithmeticError("no convergence after %s steps."%iter)
622     #
623     # return u,p
624 gross 873
625 jgs 149 # vim: expandtab shiftwidth=4:

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26