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

Contents of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1703 - (show annotations)
Thu Aug 14 05:34:25 2008 UTC (11 years, 1 month ago) by artak
File MIME type: text/x-python
File size: 112886 byte(s)
TFQMR solver is added to PASO solver. It is not parallelised yet.
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 The module provides an interface to define and solve linear partial
19 differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
20 solver capabilities in itself but hands the PDE over to
21 the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
22 The general interface is provided through the L{LinearPDE} class. The
23 L{AdvectivePDE} which is derived from the L{LinearPDE} class
24 provides an interface to PDE dominated by its advective terms. The L{Poisson},
25 L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
26 classs which are also derived form the L{LinearPDE} class should be used
27 to define of solve these sepecial PDEs.
28
29 @var __author__: name of author
30 @var __copyright__: copyrights
31 @var __license__: licence agreement
32 @var __url__: url entry point on documentation
33 @var __version__: version
34 @var __date__: date of the version
35 """
36
37 import math
38 import escript
39 import util
40 import numarray
41
42 __author__="Lutz Gross, l.gross@uq.edu.au"
43 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
44 http://www.access.edu.au
45 Primary Business: Queensland, Australia"""
46 __license__="""Licensed under the Open Software License version 3.0
47 http://www.opensource.org/licenses/osl-3.0.php"""
48 __url__="http://www.iservo.edu.au/esys"
49 __version__="$Revision$"
50 __date__="$Date$"
51
52
53 class IllegalCoefficient(ValueError):
54 """
55 raised if an illegal coefficient of the general ar particular PDE is requested.
56 """
57 pass
58
59 class IllegalCoefficientValue(ValueError):
60 """
61 raised if an incorrect value for a coefficient is used.
62 """
63 pass
64
65 class IllegalCoefficientFunctionSpace(ValueError):
66 """
67 raised if an incorrect function space for a coefficient is used.
68 """
69
70 class UndefinedPDEError(ValueError):
71 """
72 raised if a PDE is not fully defined yet.
73 """
74 pass
75
76 class PDECoefficient(object):
77 """
78 A class for describing a PDE coefficient
79
80 @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
81 @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
82 @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
83 @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
84 @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
85 @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
86 @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
87 @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
88 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
89 @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
90 @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
91 @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
92 @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
93 @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
94
95 """
96 INTERIOR=0
97 BOUNDARY=1
98 CONTACT=2
99 SOLUTION=3
100 REDUCED=4
101 BY_EQUATION=5
102 BY_SOLUTION=6
103 BY_DIM=7
104 OPERATOR=10
105 RIGHTHANDSIDE=11
106 BOTH=12
107 INTERIOR_REDUCED=13
108 BOUNDARY_REDUCED=14
109 CONTACT_REDUCED=15
110
111 def __init__(self, where, pattern, altering):
112 """
113 Initialise a PDE Coefficient type
114
115 @param where: describes where the coefficient lives
116 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
117 L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
118 @param pattern: describes the shape of the coefficient and how the shape is build for a given
119 spatial dimension and numbers of equation and solution in then PDE. For instance,
120 (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
121 is instanciated as shape (3,2,2) in case of a three equations and two solution components
122 on a 2-dimensional domain. In the case of single equation and a single solution component
123 the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
124 the example would be read as (2,).
125 @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
126 @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
127 @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
128 @param reduced: indicates if reduced
129 @type reduced: C{bool}
130 """
131 super(PDECoefficient, self).__init__()
132 self.what=where
133 self.pattern=pattern
134 self.altering=altering
135 self.resetValue()
136
137 def resetValue(self):
138 """
139 resets coefficient value to default
140 """
141 self.value=escript.Data()
142
143 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
144 """
145 defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
146
147 @param domain: domain on which the PDE uses the coefficient
148 @type domain: L{Domain<escript.Domain>}
149 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
150 @type reducedEquationOrder: C{bool}
151 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
152 @type reducedSolutionOrder: C{bool}
153 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
154 @rtype: L{FunctionSpace<escript.FunctionSpace>}
155 """
156 if self.what==self.INTERIOR:
157 return escript.Function(domain)
158 elif self.what==self.INTERIOR_REDUCED:
159 return escript.ReducedFunction(domain)
160 elif self.what==self.BOUNDARY:
161 return escript.FunctionOnBoundary(domain)
162 elif self.what==self.BOUNDARY_REDUCED:
163 return escript.ReducedFunctionOnBoundary(domain)
164 elif self.what==self.CONTACT:
165 return escript.FunctionOnContactZero(domain)
166 elif self.what==self.CONTACT_REDUCED:
167 return escript.ReducedFunctionOnContactZero(domain)
168 elif self.what==self.SOLUTION:
169 if reducedEquationOrder and reducedSolutionOrder:
170 return escript.ReducedSolution(domain)
171 else:
172 return escript.Solution(domain)
173 elif self.what==self.REDUCED:
174 return escript.ReducedSolution(domain)
175
176 def getValue(self):
177 """
178 returns the value of the coefficient
179
180 @return: value of the coefficient
181 @rtype: L{Data<escript.Data>}
182 """
183 return self.value
184
185 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
186 """
187 set the value of the coefficient to a new value
188
189 @param domain: domain on which the PDE uses the coefficient
190 @type domain: L{Domain<escript.Domain>}
191 @param numEquations: number of equations of the PDE
192 @type numEquations: C{int}
193 @param numSolutions: number of components of the PDE solution
194 @type numSolutions: C{int}
195 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
196 @type reducedEquationOrder: C{bool}
197 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
198 @type reducedSolutionOrder: C{bool}
199 @param newValue: number of components of the PDE solution
200 @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
201 @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
202 @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
203 """
204 if newValue==None:
205 newValue=escript.Data()
206 elif isinstance(newValue,escript.Data):
207 if not newValue.isEmpty():
208 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
209 try:
210 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
211 except:
212 raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
213 else:
214 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
215 if not newValue.isEmpty():
216 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
217 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
218 self.value=newValue
219
220 def isAlteringOperator(self):
221 """
222 checks if the coefficient alters the operator of the PDE
223
224 @return: True if the operator of the PDE is changed when the coefficient is changed
225 @rtype: C{bool}
226 """
227 if self.altering==self.OPERATOR or self.altering==self.BOTH:
228 return not None
229 else:
230 return None
231
232 def isAlteringRightHandSide(self):
233 """
234 checks if the coefficeint alters the right hand side of the PDE
235
236 @rtype: C{bool}
237 @return: True if the right hand side of the PDE is changed when the coefficient is changed
238 """
239 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
240 return not None
241 else:
242 return None
243
244 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
245 """
246 tries to estimate the number of equations and number of solutions if the coefficient has the given shape
247
248 @param domain: domain on which the PDE uses the coefficient
249 @type domain: L{Domain<escript.Domain>}
250 @param shape: suggested shape of the coefficient
251 @type shape: C{tuple} of C{int} values
252 @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
253 If no appropriate numbers could be identified, C{None} is returned
254 @rtype: C{tuple} of two C{int} values or C{None}
255 """
256 dim=domain.getDim()
257 if len(shape)>0:
258 num=max(shape)+1
259 else:
260 num=1
261 search=[]
262 if self.definesNumEquation() and self.definesNumSolutions():
263 for u in range(num):
264 for e in range(num):
265 search.append((e,u))
266 search.sort(self.__CompTuple2)
267 for item in search:
268 s=self.getShape(domain,item[0],item[1])
269 if len(s)==0 and len(shape)==0:
270 return (1,1)
271 else:
272 if s==shape: return item
273 elif self.definesNumEquation():
274 for e in range(num,0,-1):
275 s=self.getShape(domain,e,0)
276 if len(s)==0 and len(shape)==0:
277 return (1,None)
278 else:
279 if s==shape: return (e,None)
280
281 elif self.definesNumSolutions():
282 for u in range(num,0,-1):
283 s=self.getShape(domain,0,u)
284 if len(s)==0 and len(shape)==0:
285 return (None,1)
286 else:
287 if s==shape: return (None,u)
288 return None
289 def definesNumSolutions(self):
290 """
291 checks if the coefficient allows to estimate the number of solution components
292
293 @return: True if the coefficient allows an estimate of the number of solution components
294 @rtype: C{bool}
295 """
296 for i in self.pattern:
297 if i==self.BY_SOLUTION: return True
298 return False
299
300 def definesNumEquation(self):
301 """
302 checks if the coefficient allows to estimate the number of equations
303
304 @return: True if the coefficient allows an estimate of the number of equations
305 @rtype: C{bool}
306 """
307 for i in self.pattern:
308 if i==self.BY_EQUATION: return True
309 return False
310
311 def __CompTuple2(self,t1,t2):
312 """
313 Compare two tuples of possible number of equations and number of solutions
314
315 @param t1: The first tuple
316 @param t2: The second tuple
317
318 """
319
320 dif=t1[0]+t1[1]-(t2[0]+t2[1])
321 if dif<0: return 1
322 elif dif>0: return -1
323 else: return 0
324
325 def getShape(self,domain,numEquations=1,numSolutions=1):
326 """
327 builds the required shape of the coefficient
328
329 @param domain: domain on which the PDE uses the coefficient
330 @type domain: L{Domain<escript.Domain>}
331 @param numEquations: number of equations of the PDE
332 @type numEquations: C{int}
333 @param numSolutions: number of components of the PDE solution
334 @type numSolutions: C{int}
335 @return: shape of the coefficient
336 @rtype: C{tuple} of C{int} values
337 """
338 dim=domain.getDim()
339 s=()
340 for i in self.pattern:
341 if i==self.BY_EQUATION:
342 if numEquations>1: s=s+(numEquations,)
343 elif i==self.BY_SOLUTION:
344 if numSolutions>1: s=s+(numSolutions,)
345 else:
346 s=s+(dim,)
347 return s
348
349 class LinearPDE(object):
350 """
351 This class is used to define a general linear, steady, second order PDE
352 for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
353
354 For a single PDE with a solution with a single component the linear PDE is defined in the following form:
355
356 M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
357
358
359 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
360 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
361 The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
362 L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
363 L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
364 such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B_reduced}, M{C_reduced}, M{X_reduced}
365 M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} and M{Y_reduced} are scalar.
366
367 The following natural boundary conditions are considered:
368
369 M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
370
371 where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
372
373
374 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
375
376 M{u=r} where M{q>0}
377
378 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
379 The constraints override any other condition set by the PDE or the boundary condition.
380
381 The PDE is symmetrical if
382
383 M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j]}
384
385 For a system of PDEs and a solution with several components the PDE has the form
386
387 M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
388
389 M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.
390 The natural boundary conditions take the form:
391
392 M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
393
394
395 The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
396
397 Constraints take the form
398
399 M{u[i]=r[i]} where M{q[i]>0}
400
401 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
402
403 The system of PDEs is symmetrical if
404
405 - M{A[i,j,k,l]=A[k,l,i,j]}
406 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
407 - M{B[i,j,k]=C[k,i,j]}
408 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
409 - M{D[i,k]=D[i,k]}
410 - M{D_reduced[i,k]=D_reduced[i,k]}
411 - M{d[i,k]=d[k,i]}
412 - M{d_reduced[i,k]=d_reduced[k,i]}
413
414 L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
415 discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
416 defined as
417
418 M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
419
420 For the case of single solution component and single PDE M{J} is defined
421
422 M{J_{j}=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
423
424 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
425 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
426 the contact condition takes the form
427
428 M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
429
430 where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
431 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
432 L{jump<util.jump>}.
433 The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
434 The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
435 In case of a single PDE and a single component solution the contact condition takes the form
436
437 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
438
439 In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
440
441 @cvar DEFAULT: The default method used to solve the system of linear equations
442 @cvar DIRECT: The direct solver based on LDU factorization
443 @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
444 @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
445 @cvar CR: The conjugate residual method
446 @cvar CGS: The conjugate gardient square method
447 @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
448 @cvar TFQMR: Transport Free Quasi Minimal Residual method.
449 @cvar SSOR: The symmetric overrealaxtion method
450 @cvar ILU0: The incomplete LU factorization preconditioner with no fill in
451 @cvar ILUT: The incomplete LU factorization preconditioner with will in
452 @cvar JACOBI: The Jacobi preconditioner
453 @cvar GMRES: The Gram-Schmidt minimum residual method
454 @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
455 @cvar LUMPING: Matrix lumping.
456 @cvar NO_REORDERING: No matrix reordering allowed
457 @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
458 @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
459 @cvar PASO: PASO solver package
460 @cvar SCSL: SGI SCSL solver library
461 @cvar MKL: Intel's MKL solver library
462 @cvar UMFPACK: the UMFPACK library
463 @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
464 @cvar ITERATIVE: The default iterative solver
465 @cvar AMG: algebraic multi grid
466 @cvar RILU: recursive ILU
467
468 """
469 DEFAULT= 0
470 DIRECT= 1
471 CHOLEVSKY= 2
472 PCG= 3
473 CR= 4
474 CGS= 5
475 BICGSTAB= 6
476 SSOR= 7
477 ILU0= 8
478 ILUT= 9
479 JACOBI= 10
480 GMRES= 11
481 PRES20= 12
482 LUMPING= 13
483 NO_REORDERING= 17
484 MINIMUM_FILL_IN= 18
485 NESTED_DISSECTION= 19
486 SCSL= 14
487 MKL= 15
488 UMFPACK= 16
489 ITERATIVE= 20
490 PASO= 21
491 AMG= 22
492 RILU = 23
493 TRILINOS = 24
494 NONLINEAR_GMRES = 25
495 TFQMR = 26
496
497 SMALL_TOLERANCE=1.e-13
498 __PACKAGE_KEY="package"
499 __METHOD_KEY="method"
500 __SYMMETRY_KEY="symmetric"
501 __TOLERANCE_KEY="tolerance"
502 __PRECONDITIONER_KEY="preconditioner"
503
504
505 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
506 """
507 initializes a new linear PDE
508
509 @param domain: domain of the PDE
510 @type domain: L{Domain<escript.Domain>}
511 @param numEquations: number of equations. If numEquations==None the number of equations
512 is exracted from the PDE coefficients.
513 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
514 is exracted from the PDE coefficients.
515 @param debug: if True debug informations are printed.
516
517 """
518 super(LinearPDE, self).__init__()
519 #
520 # the coefficients of the general PDE:
521 #
522 self.__COEFFICIENTS_OF_GENEARL_PDE={
523 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
524 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
525 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
526 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
527 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
528 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
529 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
530 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
531 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
532 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
533 "A_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
534 "B_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
535 "C_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
536 "D_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
537 "X_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
538 "Y_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
539 "d_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
540 "y_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
541 "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
542 "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
543 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
544 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
545
546 # COEFFICIENTS can be overwritten by subclasses:
547 self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
548 self.__altered_coefficients=False
549 # initialize attributes
550 self.__debug=debug
551 self.__domain=domain
552 self.__numEquations=numEquations
553 self.__numSolutions=numSolutions
554 self.__resetSystem()
555
556 # set some default values:
557 self.__reduce_equation_order=False
558 self.__reduce_solution_order=False
559 self.__tolerance=1.e-8
560 self.__solver_method=self.DEFAULT
561 self.__solver_package=self.DEFAULT
562 self.__preconditioner=self.DEFAULT
563 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
564 self.__sym=False
565
566 self.resetCoefficients()
567 self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
568 # =============================================================================
569 # general stuff:
570 # =============================================================================
571 def __str__(self):
572 """
573 returns string representation of the PDE
574
575 @return: a simple representation of the PDE
576 @rtype: C{str}
577 """
578 return "<LinearPDE %d>"%id(self)
579 # =============================================================================
580 # debug :
581 # =============================================================================
582 def setDebugOn(self):
583 """
584 switches on debugging
585 """
586 self.__debug=not None
587
588 def setDebugOff(self):
589 """
590 switches off debugging
591 """
592 self.__debug=None
593
594 def trace(self,text):
595 """
596 print the text message if debugging is swiched on.
597 @param text: message
598 @type text: C{string}
599 """
600 if self.__debug: print "%s: %s"%(str(self),text)
601
602 # =============================================================================
603 # some service functions:
604 # =============================================================================
605 def getDomain(self):
606 """
607 returns the domain of the PDE
608
609 @return: the domain of the PDE
610 @rtype: L{Domain<escript.Domain>}
611 """
612 return self.__domain
613
614 def getDim(self):
615 """
616 returns the spatial dimension of the PDE
617
618 @return: the spatial dimension of the PDE domain
619 @rtype: C{int}
620 """
621 return self.getDomain().getDim()
622
623 def getNumEquations(self):
624 """
625 returns the number of equations
626
627 @return: the number of equations
628 @rtype: C{int}
629 @raise UndefinedPDEError: if the number of equations is not be specified yet.
630 """
631 if self.__numEquations==None:
632 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
633 else:
634 return self.__numEquations
635
636 def getNumSolutions(self):
637 """
638 returns the number of unknowns
639
640 @return: the number of unknowns
641 @rtype: C{int}
642 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
643 """
644 if self.__numSolutions==None:
645 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
646 else:
647 return self.__numSolutions
648
649 def reduceEquationOrder(self):
650 """
651 return status for order reduction for equation
652
653 @return: return True is reduced interpolation order is used for the represenation of the equation
654 @rtype: L{bool}
655 """
656 return self.__reduce_equation_order
657
658 def reduceSolutionOrder(self):
659 """
660 return status for order reduction for the solution
661
662 @return: return True is reduced interpolation order is used for the represenation of the solution
663 @rtype: L{bool}
664 """
665 return self.__reduce_solution_order
666
667 def getFunctionSpaceForEquation(self):
668 """
669 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
670
671 @return: representation space of equation
672 @rtype: L{FunctionSpace<escript.FunctionSpace>}
673 """
674 if self.reduceEquationOrder():
675 return escript.ReducedSolution(self.getDomain())
676 else:
677 return escript.Solution(self.getDomain())
678
679 def getFunctionSpaceForSolution(self):
680 """
681 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
682
683 @return: representation space of solution
684 @rtype: L{FunctionSpace<escript.FunctionSpace>}
685 """
686 if self.reduceSolutionOrder():
687 return escript.ReducedSolution(self.getDomain())
688 else:
689 return escript.Solution(self.getDomain())
690
691
692 def getOperator(self):
693 """
694 provides access to the operator of the PDE
695
696 @return: the operator of the PDE
697 @rtype: L{Operator<escript.Operator>}
698 """
699 m=self.getSystem()[0]
700 if self.isUsingLumping():
701 return self.copyConstraint(1./m)
702 else:
703 return m
704
705 def getRightHandSide(self):
706 """
707 provides access to the right hand side of the PDE
708 @return: the right hand side of the PDE
709 @rtype: L{Data<escript.Data>}
710 """
711 r=self.getSystem()[1]
712 if self.isUsingLumping():
713 return self.copyConstraint(r)
714 else:
715 return r
716
717 def applyOperator(self,u=None):
718 """
719 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
720
721 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
722 the current solution is used.
723 @type u: L{Data<escript.Data>} or None
724 @return: image of u
725 @rtype: L{Data<escript.Data>}
726 """
727 if u==None:
728 return self.getOperator()*self.getSolution()
729 else:
730 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
731
732 def getResidual(self,u=None):
733 """
734 return the residual of u or the current solution if u is not present.
735
736 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
737 the current solution is used.
738 @type u: L{Data<escript.Data>} or None
739 @return: residual of u
740 @rtype: L{Data<escript.Data>}
741 """
742 return self.applyOperator(u)-self.getRightHandSide()
743
744 def checkSymmetry(self,verbose=True):
745 """
746 test the PDE for symmetry.
747
748 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
749 @type verbose: C{bool}
750 @return: True if the PDE is symmetric.
751 @rtype: L{Data<escript.Data>}
752 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
753 """
754 verbose=verbose or self.__debug
755 out=True
756 if self.getNumSolutions()!=self.getNumEquations():
757 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
758 out=False
759 else:
760 A=self.getCoefficientOfGeneralPDE("A")
761 if not A.isEmpty():
762 tol=util.Lsup(A)*self.SMALL_TOLERANCE
763 if self.getNumSolutions()>1:
764 for i in range(self.getNumEquations()):
765 for j in range(self.getDim()):
766 for k in range(self.getNumSolutions()):
767 for l in range(self.getDim()):
768 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
769 if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
770 out=False
771 else:
772 for j in range(self.getDim()):
773 for l in range(self.getDim()):
774 if util.Lsup(A[j,l]-A[l,j])>tol:
775 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
776 out=False
777 B=self.getCoefficientOfGeneralPDE("B")
778 C=self.getCoefficientOfGeneralPDE("C")
779 if B.isEmpty() and not C.isEmpty():
780 if verbose: print "non-symmetric PDE because B is not present but C is"
781 out=False
782 elif not B.isEmpty() and C.isEmpty():
783 if verbose: print "non-symmetric PDE because C is not present but B is"
784 out=False
785 elif not B.isEmpty() and not C.isEmpty():
786 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
787 if self.getNumSolutions()>1:
788 for i in range(self.getNumEquations()):
789 for j in range(self.getDim()):
790 for k in range(self.getNumSolutions()):
791 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
792 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
793 out=False
794 else:
795 for j in range(self.getDim()):
796 if util.Lsup(B[j]-C[j])>tol:
797 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
798 out=False
799 if self.getNumSolutions()>1:
800 D=self.getCoefficientOfGeneralPDE("D")
801 if not D.isEmpty():
802 tol=util.Lsup(D)*self.SMALL_TOLERANCE
803 for i in range(self.getNumEquations()):
804 for k in range(self.getNumSolutions()):
805 if util.Lsup(D[i,k]-D[k,i])>tol:
806 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
807 out=False
808 d=self.getCoefficientOfGeneralPDE("d")
809 if not d.isEmpty():
810 tol=util.Lsup(d)*self.SMALL_TOLERANCE
811 for i in range(self.getNumEquations()):
812 for k in range(self.getNumSolutions()):
813 if util.Lsup(d[i,k]-d[k,i])>tol:
814 if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
815 out=False
816 d_contact=self.getCoefficientOfGeneralPDE("d_contact")
817 if not d_contact.isEmpty():
818 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
819 for i in range(self.getNumEquations()):
820 for k in range(self.getNumSolutions()):
821 if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
822 if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
823 out=False
824 # and now the reduced coefficients
825 A_reduced=self.getCoefficientOfGeneralPDE("A_reduced")
826 if not A_reduced.isEmpty():
827 tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE
828 if self.getNumSolutions()>1:
829 for i in range(self.getNumEquations()):
830 for j in range(self.getDim()):
831 for k in range(self.getNumSolutions()):
832 for l in range(self.getDim()):
833 if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol:
834 if verbose: print "non-symmetric PDE because A_reduced[%d,%d,%d,%d]!=A_reduced[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
835 out=False
836 else:
837 for j in range(self.getDim()):
838 for l in range(self.getDim()):
839 if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol:
840 if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j)
841 out=False
842 B_reduced=self.getCoefficientOfGeneralPDE("B_reduced")
843 C_reduced=self.getCoefficientOfGeneralPDE("C_reduced")
844 if B_reduced.isEmpty() and not C_reduced.isEmpty():
845 if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is"
846 out=False
847 elif not B_reduced.isEmpty() and C_reduced.isEmpty():
848 if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is"
849 out=False
850 elif not B_reduced.isEmpty() and not C_reduced.isEmpty():
851 tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2.
852 if self.getNumSolutions()>1:
853 for i in range(self.getNumEquations()):
854 for j in range(self.getDim()):
855 for k in range(self.getNumSolutions()):
856 if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol:
857 if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j)
858 out=False
859 else:
860 for j in range(self.getDim()):
861 if util.Lsup(B_reduced[j]-C_reduced[j])>tol:
862 if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j)
863 out=False
864 if self.getNumSolutions()>1:
865 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
866 if not D_reduced.isEmpty():
867 tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE
868 for i in range(self.getNumEquations()):
869 for k in range(self.getNumSolutions()):
870 if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol:
871 if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i)
872 out=False
873 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
874 if not d_reduced.isEmpty():
875 tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE
876 for i in range(self.getNumEquations()):
877 for k in range(self.getNumSolutions()):
878 if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol:
879 if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i)
880 out=False
881 d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
882 if not d_contact_reduced.isEmpty():
883 tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE
884 for i in range(self.getNumEquations()):
885 for k in range(self.getNumSolutions()):
886 if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol:
887 if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i)
888 out=False
889 return out
890
891 def getSolution(self,**options):
892 """
893 returns the solution of the PDE. If the solution is not valid the PDE is solved.
894
895 @return: the solution
896 @rtype: L{Data<escript.Data>}
897 @param options: solver options
898 @keyword verbose: True to get some information during PDE solution
899 @type verbose: C{bool}
900 @keyword reordering: reordering scheme to be used during elimination. Allowed values are
901 L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
902 @keyword iter_max: maximum number of iteration steps allowed.
903 @keyword drop_tolerance: threshold for drupping in L{ILUT}
904 @keyword drop_storage: maximum of allowed memory in L{ILUT}
905 @keyword truncation: maximum number of residuals in L{GMRES}
906 @keyword restart: restart cycle length in L{GMRES}
907 """
908 if not self.__solution_isValid:
909 mat,f=self.getSystem()
910 if self.isUsingLumping():
911 self.__solution=self.copyConstraint(f*mat)
912 else:
913 options[self.__TOLERANCE_KEY]=self.getTolerance()
914 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
915 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
916 options[self.__PACKAGE_KEY]=self.getSolverPackage()
917 options[self.__SYMMETRY_KEY]=self.isSymmetric()
918 self.trace("PDE is resolved.")
919 self.trace("solver options: %s"%str(options))
920 self.__solution=mat.solve(f,options)
921 self.__solution_isValid=True
922 return self.__solution
923
924 def getFlux(self,u=None):
925 """
926 returns the flux M{J} for a given M{u}
927
928 M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}
929
930 or
931
932 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
933
934 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
935 @type u: L{Data<escript.Data>} or None
936 @return: flux
937 @rtype: L{Data<escript.Data>}
938 """
939 if u==None: u=self.getSolution()
940 return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \
941 +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \
942 -util.self.getCoefficientOfGeneralPDE("X") \
943 +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
944 +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \
945 -util.self.getCoefficientOfGeneralPDE("X_reduced")
946 # =============================================================================
947 # solver settings:
948 # =============================================================================
949 def setSolverMethod(self,solver=None,preconditioner=None):
950 """
951 sets a new solver
952
953 @param solver: sets a new solver method.
954 @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{PRES20}, L{LUMPING}, L{AMG}
955 @param preconditioner: sets a new solver method.
956 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
957 """
958 if solver==None: solver=self.__solver_method
959 if preconditioner==None: preconditioner=self.__preconditioner
960 if solver==None: solver=self.DEFAULT
961 if preconditioner==None: preconditioner=self.DEFAULT
962 if not (solver,preconditioner)==self.getSolverMethod():
963 self.__solver_method=solver
964 self.__preconditioner=preconditioner
965 self.__checkMatrixType()
966 self.trace("New solver is %s"%self.getSolverMethodName())
967
968 def getSolverMethodName(self):
969 """
970 returns the name of the solver currently used
971
972 @return: the name of the solver currently used.
973 @rtype: C{string}
974 """
975
976 m=self.getSolverMethod()
977 p=self.getSolverPackage()
978 method=""
979 if m[0]==self.DEFAULT: method="DEFAULT"
980 elif m[0]==self.DIRECT: method= "DIRECT"
981 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
982 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
983 elif m[0]==self.PCG: method= "PCG"
984 elif m[0]==self.PCG: method= "TFQMR"
985 elif m[0]==self.CR: method= "CR"
986 elif m[0]==self.CGS: method= "CGS"
987 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
988 elif m[0]==self.SSOR: method= "SSOR"
989 elif m[0]==self.GMRES: method= "GMRES"
990 elif m[0]==self.PRES20: method= "PRES20"
991 elif m[0]==self.LUMPING: method= "LUMPING"
992 elif m[0]==self.AMG: method= "AMG"
993 if m[1]==self.DEFAULT: method+="+DEFAULT"
994 elif m[1]==self.JACOBI: method+= "+JACOBI"
995 elif m[1]==self.ILU0: method+= "+ILU0"
996 elif m[1]==self.ILUT: method+= "+ILUT"
997 elif m[1]==self.SSOR: method+= "+SSOR"
998 elif m[1]==self.AMG: method+= "+AMG"
999 elif m[1]==self.RILU: method+= "+RILU"
1000 if p==self.DEFAULT: package="DEFAULT"
1001 elif p==self.PASO: package= "PASO"
1002 elif p==self.MKL: package= "MKL"
1003 elif p==self.SCSL: package= "SCSL"
1004 elif p==self.UMFPACK: package= "UMFPACK"
1005 elif p==self.TRILINOS: package= "TRILINOS"
1006 else : method="unknown"
1007 return "%s solver of %s package"%(method,package)
1008
1009
1010 def getSolverMethod(self):
1011 """
1012 returns the solver method
1013
1014 @return: the solver method currently be used.
1015 @rtype: C{int}
1016 """
1017 return self.__solver_method,self.__preconditioner
1018
1019 def setSolverPackage(self,package=None):
1020 """
1021 sets a new solver package
1022
1023 @param package: sets a new solver method.
1024 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1025 """
1026 if package==None: package=self.DEFAULT
1027 if not package==self.getSolverPackage():
1028 self.__solver_package=package
1029 self.__checkMatrixType()
1030 self.trace("New solver is %s"%self.getSolverMethodName())
1031
1032 def getSolverPackage(self):
1033 """
1034 returns the package of the solver
1035
1036 @return: the solver package currently being used.
1037 @rtype: C{int}
1038 """
1039 return self.__solver_package
1040
1041 def isUsingLumping(self):
1042 """
1043 checks if matrix lumping is used a solver method
1044
1045 @return: True is lumping is currently used a solver method.
1046 @rtype: C{bool}
1047 """
1048 return self.getSolverMethod()[0]==self.LUMPING
1049
1050 def setTolerance(self,tol=1.e-8):
1051 """
1052 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
1053
1054 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
1055
1056 defines the stopping criterion.
1057
1058 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
1059 the system will be resolved.
1060 @type tol: positive C{float}
1061 @raise ValueError: if tolerance is not positive.
1062 """
1063 if not tol>0:
1064 raise ValueError,"Tolerance as to be positive"
1065 if tol<self.getTolerance(): self.__invalidateSolution()
1066 self.trace("New tolerance %e"%tol)
1067 self.__tolerance=tol
1068 return
1069
1070 def getTolerance(self):
1071 """
1072 returns the tolerance set for the solution
1073
1074 @return: tolerance currently used.
1075 @rtype: C{float}
1076 """
1077 return self.__tolerance
1078
1079 # =============================================================================
1080 # symmetry flag:
1081 # =============================================================================
1082 def isSymmetric(self):
1083 """
1084 checks if symmetry is indicated.
1085
1086 @return: True is a symmetric PDE is indicated, otherwise False is returned
1087 @rtype: C{bool}
1088 """
1089 return self.__sym
1090
1091 def setSymmetryOn(self):
1092 """
1093 sets the symmetry flag.
1094 """
1095 if not self.isSymmetric():
1096 self.trace("PDE is set to be symmetric")
1097 self.__sym=True
1098 self.__checkMatrixType()
1099
1100 def setSymmetryOff(self):
1101 """
1102 removes the symmetry flag.
1103 """
1104 if self.isSymmetric():
1105 self.trace("PDE is set to be unsymmetric")
1106 self.__sym=False
1107 self.__checkMatrixType()
1108
1109 def setSymmetryTo(self,flag=False):
1110 """
1111 sets the symmetry flag to flag
1112
1113 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
1114 @type flag: C{bool}
1115 """
1116 if flag:
1117 self.setSymmetryOn()
1118 else:
1119 self.setSymmetryOff()
1120
1121 # =============================================================================
1122 # function space handling for the equation as well as the solution
1123 # =============================================================================
1124 def setReducedOrderOn(self):
1125 """
1126 switches on reduced order for solution and equation representation
1127
1128 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1129 """
1130 self.setReducedOrderForSolutionOn()
1131 self.setReducedOrderForEquationOn()
1132
1133 def setReducedOrderOff(self):
1134 """
1135 switches off reduced order for solution and equation representation
1136
1137 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1138 """
1139 self.setReducedOrderForSolutionOff()
1140 self.setReducedOrderForEquationOff()
1141
1142 def setReducedOrderTo(self,flag=False):
1143 """
1144 sets order reduction for both solution and equation representation according to flag.
1145 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
1146 if flag is not present order reduction is switched off
1147 @type flag: C{bool}
1148 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1149 """
1150 self.setReducedOrderForSolutionTo(flag)
1151 self.setReducedOrderForEquationTo(flag)
1152
1153
1154 def setReducedOrderForSolutionOn(self):
1155 """
1156 switches on reduced order for solution representation
1157
1158 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1159 """
1160 if not self.__reduce_solution_order:
1161 if self.__altered_coefficients:
1162 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1163 self.trace("Reduced order is used to solution representation.")
1164 self.__reduce_solution_order=True
1165 self.__resetSystem()
1166
1167 def setReducedOrderForSolutionOff(self):
1168 """
1169 switches off reduced order for solution representation
1170
1171 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1172 """
1173 if self.__reduce_solution_order:
1174 if self.__altered_coefficients:
1175 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1176 self.trace("Full order is used to interpolate solution.")
1177 self.__reduce_solution_order=False
1178 self.__resetSystem()
1179
1180 def setReducedOrderForSolutionTo(self,flag=False):
1181 """
1182 sets order for test functions according to flag
1183
1184 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1185 if flag is not present order reduction is switched off
1186 @type flag: C{bool}
1187 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1188 """
1189 if flag:
1190 self.setReducedOrderForSolutionOn()
1191 else:
1192 self.setReducedOrderForSolutionOff()
1193
1194 def setReducedOrderForEquationOn(self):
1195 """
1196 switches on reduced order for equation representation
1197
1198 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1199 """
1200 if not self.__reduce_equation_order:
1201 if self.__altered_coefficients:
1202 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1203 self.trace("Reduced order is used for test functions.")
1204 self.__reduce_equation_order=True
1205 self.__resetSystem()
1206
1207 def setReducedOrderForEquationOff(self):
1208 """
1209 switches off reduced order for equation representation
1210
1211 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1212 """
1213 if self.__reduce_equation_order:
1214 if self.__altered_coefficients:
1215 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1216 self.trace("Full order is used for test functions.")
1217 self.__reduce_equation_order=False
1218 self.__resetSystem()
1219
1220 def setReducedOrderForEquationTo(self,flag=False):
1221 """
1222 sets order for test functions according to flag
1223
1224 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1225 if flag is not present order reduction is switched off
1226 @type flag: C{bool}
1227 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1228 """
1229 if flag:
1230 self.setReducedOrderForEquationOn()
1231 else:
1232 self.setReducedOrderForEquationOff()
1233
1234 # =============================================================================
1235 # private method:
1236 # =============================================================================
1237 def __checkMatrixType(self):
1238 """
1239 reassess the matrix type and, if a new matrix is needed, resets the system.
1240 """
1241 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1242 if not new_matrix_type==self.__matrix_type:
1243 self.trace("Matrix type is now %d."%new_matrix_type)
1244 self.__matrix_type=new_matrix_type
1245 self.__resetSystem()
1246 #
1247 # rebuild switches :
1248 #
1249 def __invalidateSolution(self):
1250 """
1251 indicates the PDE has to be resolved if the solution is requested
1252 """
1253 if self.__solution_isValid: self.trace("PDE has to be resolved.")
1254 self.__solution_isValid=False
1255
1256 def __invalidateOperator(self):
1257 """
1258 indicates the operator has to be rebuilt next time it is used
1259 """
1260 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1261 self.__invalidateSolution()
1262 self.__operator_is_Valid=False
1263
1264 def __invalidateRightHandSide(self):
1265 """
1266 indicates the right hand side has to be rebuild next time it is used
1267 """
1268 if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1269 self.__invalidateSolution()
1270 self.__righthandside_isValid=False
1271
1272 def __invalidateSystem(self):
1273 """
1274 annonced that everthing has to be rebuild:
1275 """
1276 if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1277 self.__invalidateSolution()
1278 self.__invalidateOperator()
1279 self.__invalidateRightHandSide()
1280
1281 def __resetSystem(self):
1282 """
1283 annonced that everthing has to be rebuild:
1284 """
1285 self.trace("New System is built from scratch.")
1286 self.__operator=escript.Operator()
1287 self.__operator_is_Valid=False
1288 self.__righthandside=escript.Data()
1289 self.__righthandside_isValid=False
1290 self.__solution=escript.Data()
1291 self.__solution_isValid=False
1292 #
1293 # system initialization:
1294 #
1295 def __getNewOperator(self):
1296 """
1297 returns an instance of a new operator
1298 """
1299 self.trace("New operator is allocated.")
1300 return self.getDomain().newOperator( \
1301 self.getNumEquations(), \
1302 self.getFunctionSpaceForEquation(), \
1303 self.getNumSolutions(), \
1304 self.getFunctionSpaceForSolution(), \
1305 self.__matrix_type)
1306
1307 def __getNewRightHandSide(self):
1308 """
1309 returns an instance of a new right hand side
1310 """
1311 self.trace("New right hand side is allocated.")
1312 if self.getNumEquations()>1:
1313 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1314 else:
1315 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1316
1317 def __getNewSolution(self):
1318 """
1319 returns an instance of a new solution
1320 """
1321 self.trace("New solution is allocated.")
1322 if self.getNumSolutions()>1:
1323 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1324 else:
1325 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1326
1327 def __makeFreshSolution(self):
1328 """
1329 makes sure that the solution is instantiated and returns it initialized by zeros
1330 """
1331 if self.__solution.isEmpty():
1332 self.__solution=self.__getNewSolution()
1333 else:
1334 self.__solution*=0
1335 self.trace("Solution is reset to zero.")
1336 return self.__solution
1337
1338 def __makeFreshRightHandSide(self):
1339 """
1340 makes sure that the right hand side is instantiated and returns it initialized by zeros
1341 """
1342 if self.__righthandside.isEmpty():
1343 self.__righthandside=self.__getNewRightHandSide()
1344 else:
1345 self.__righthandside.setToZero()
1346 self.trace("Right hand side is reset to zero.")
1347 return self.__righthandside
1348
1349 def __makeFreshOperator(self):
1350 """
1351 makes sure that the operator is instantiated and returns it initialized by zeros
1352 """
1353 if self.__operator.isEmpty():
1354 self.__operator=self.__getNewOperator()
1355 else:
1356 self.__operator.resetValues()
1357 self.trace("Operator reset to zero")
1358 return self.__operator
1359
1360 def __applyConstraint(self):
1361 """
1362 applies the constraints defined by q and r to the system
1363 """
1364 if not self.isUsingLumping():
1365 q=self.getCoefficientOfGeneralPDE("q")
1366 r=self.getCoefficientOfGeneralPDE("r")
1367 if not q.isEmpty() and not self.__operator.isEmpty():
1368 # q is the row and column mask to indicate where constraints are set:
1369 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1370 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1371 u=self.__getNewSolution()
1372 if r.isEmpty():
1373 r_s=self.__getNewSolution()
1374 else:
1375 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1376 u.copyWithMask(r_s,col_q)
1377 if not self.__righthandside.isEmpty():
1378 self.__righthandside-=self.__operator*u
1379 self.__righthandside=self.copyConstraint(self.__righthandside)
1380 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1381 # =============================================================================
1382 # function giving access to coefficients of the general PDE:
1383 # =============================================================================
1384 def getCoefficientOfGeneralPDE(self,name):
1385 """
1386 return the value of the coefficient name of the general PDE.
1387
1388 @note: This method is called by the assembling routine it can be overwritten
1389 to map coefficients of a particular PDE to the general PDE.
1390 @param name: name of the coefficient requested.
1391 @type name: C{string}
1392 @return: the value of the coefficient name
1393 @rtype: L{Data<escript.Data>}
1394 @raise IllegalCoefficient: if name is not one of coefficients
1395 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1396 M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1397 """
1398 if self.hasCoefficientOfGeneralPDE(name):
1399 return self.getCoefficient(name)
1400 else:
1401 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1402
1403 def hasCoefficientOfGeneralPDE(self,name):
1404 """
1405 checks if name is a the name of a coefficient of the general PDE.
1406
1407 @param name: name of the coefficient enquired.
1408 @type name: C{string}
1409 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1410 @rtype: C{bool}
1411
1412 """
1413 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1414
1415 def createCoefficientOfGeneralPDE(self,name):
1416 """
1417 returns a new instance of a coefficient for coefficient name of the general PDE
1418
1419 @param name: name of the coefficient requested.
1420 @type name: C{string}
1421 @return: a coefficient name initialized to 0.
1422 @rtype: L{Data<escript.Data>}
1423 @raise IllegalCoefficient: if name is not one of coefficients
1424 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1425 M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1426 """
1427 if self.hasCoefficientOfGeneralPDE(name):
1428 return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1429 else:
1430 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1431
1432 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1433 """
1434 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1435
1436 @param name: name of the coefficient enquired.
1437 @type name: C{string}
1438 @return: the function space to be used for coefficient name
1439 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1440 @raise IllegalCoefficient: if name is not one of coefficients
1441 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1442 M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1443 """
1444 if self.hasCoefficientOfGeneralPDE(name):
1445 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1446 else:
1447 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1448
1449 def getShapeOfCoefficientOfGeneralPDE(self,name):
1450 """
1451 return the shape of the coefficient name of the general PDE
1452
1453 @param name: name of the coefficient enquired.
1454 @type name: C{string}
1455 @return: the shape of the coefficient name
1456 @rtype: C{tuple} of C{int}
1457 @raise IllegalCoefficient: if name is not one of coefficients
1458 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1459 M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced}, M{Y_reduced}, M{d_reduced}, M{y_reduced}, M{d_contact_reduced}, M{y_contact_reduced}, M{r} or M{q}.
1460 """
1461 if self.hasCoefficientOfGeneralPDE(name):
1462 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1463 else:
1464 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1465
1466 # =============================================================================
1467 # functions giving access to coefficients of a particular PDE implementation:
1468 # =============================================================================
1469 def getCoefficient(self,name):
1470 """
1471 returns the value of the coefficient name
1472
1473 @param name: name of the coefficient requested.
1474 @type name: C{string}
1475 @return: the value of the coefficient name
1476 @rtype: L{Data<escript.Data>}
1477 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1478 """
1479 if self.hasCoefficient(name):
1480 return self.COEFFICIENTS[name].getValue()
1481 else:
1482 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1483
1484 def hasCoefficient(self,name):
1485 """
1486 return True if name is the name of a coefficient
1487
1488 @param name: name of the coefficient enquired.
1489 @type name: C{string}
1490 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1491 @rtype: C{bool}
1492 """
1493 return self.COEFFICIENTS.has_key(name)
1494
1495 def createCoefficient(self, name):
1496 """
1497 create a L{Data<escript.Data>} object corresponding to coefficient name
1498
1499 @return: a coefficient name initialized to 0.
1500 @rtype: L{Data<escript.Data>}
1501 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1502 """
1503 if self.hasCoefficient(name):
1504 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1505 else:
1506 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1507
1508 def getFunctionSpaceForCoefficient(self,name):
1509 """
1510 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1511
1512 @param name: name of the coefficient enquired.
1513 @type name: C{string}
1514 @return: the function space to be used for coefficient name
1515 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1516 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1517 """
1518 if self.hasCoefficient(name):
1519 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1520 else:
1521 raise ValueError,"unknown coefficient %s requested"%name
1522 def getShapeOfCoefficient(self,name):
1523 """
1524 return the shape of the coefficient name
1525
1526 @param name: name of the coefficient enquired.
1527 @type name: C{string}
1528 @return: the shape of the coefficient name
1529 @rtype: C{tuple} of C{int}
1530 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1531 """
1532 if self.hasCoefficient(name):
1533 return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1534 else:
1535 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1536
1537 def resetCoefficients(self):
1538 """
1539 resets all coefficients to there default values.
1540 """
1541 for i in self.COEFFICIENTS.iterkeys():
1542 self.COEFFICIENTS[i].resetValue()
1543
1544 def alteredCoefficient(self,name):
1545 """
1546 announce that coefficient name has been changed
1547
1548 @param name: name of the coefficient enquired.
1549 @type name: C{string}
1550 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1551 @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.
1552 """
1553 if self.hasCoefficient(name):
1554 self.trace("Coefficient %s has been altered."%name)
1555 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1556 if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1557 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1558 else:
1559 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1560
1561 def copyConstraint(self,u):
1562 """
1563 copies the constraint into u and returns u.
1564
1565 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1566 @type u: L{Data<escript.Data>}
1567 @return: the input u modified by the constraints.
1568 @rtype: L{Data<escript.Data>}
1569 @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1570 """
1571 q=self.getCoefficientOfGeneralPDE("q")
1572 r=self.getCoefficientOfGeneralPDE("r")
1573 if not q.isEmpty():
1574 if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1575 if r.isEmpty():
1576 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1577 else:
1578 r=escript.Data(r,u.getFunctionSpace())
1579 u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1580 return u
1581
1582 def setValue(self,**coefficients):
1583 """
1584 sets new values to coefficients
1585
1586 @param coefficients: new values assigned to coefficients
1587 @keyword A: value for coefficient A.
1588 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1589 @keyword A_reduced: value for coefficient A_reduced.
1590 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1591 @keyword B: value for coefficient B
1592 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1593 @keyword B_reduced: value for coefficient B_reduced
1594 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1595 @keyword C: value for coefficient C
1596 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1597 @keyword C_reduced: value for coefficient C_reduced
1598 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1599 @keyword D: value for coefficient D
1600 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1601 @keyword D_reduced: value for coefficient D_reduced
1602 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1603 @keyword X: value for coefficient X
1604 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1605 @keyword X_reduced: value for coefficient X_reduced
1606 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1607 @keyword Y: value for coefficient Y
1608 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1609 @keyword Y_reduced: value for coefficient Y_reduced
1610 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1611 @keyword d: value for coefficient d
1612 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1613 @keyword d_reduced: value for coefficient d_reduced
1614 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1615 @keyword y: value for coefficient y
1616 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1617 @keyword d_contact: value for coefficient d_contact
1618 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1619 @keyword d_contact_reduced: value for coefficient d_contact_reduced
1620 @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
1621 @keyword y_contact: value for coefficient y_contact
1622 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1623 @keyword y_contact_reduced: value for coefficient y_contact_reduced
1624 @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
1625 @keyword r: values prescribed to the solution at the locations of constraints
1626 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1627 depending of reduced order is used for the solution.
1628 @keyword q: mask for location of constraints
1629 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1630 depending of reduced order is used for the representation of the equation.
1631 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1632 """
1633 # check if the coefficients are legal:
1634 for i in coefficients.iterkeys():
1635 if not self.hasCoefficient(i):
1636 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1637 # if the number of unknowns or equations is still unknown we try to estimate them:
1638 if self.__numEquations==None or self.__numSolutions==None:
1639 for i,d in coefficients.iteritems():
1640 if hasattr(d,"shape"):
1641 s=d.shape
1642 elif hasattr(d,"getShape"):
1643 s=d.getShape()
1644 else:
1645 s=numarray.array(d).shape
1646 if s!=None:
1647 # get number of equations and number of unknowns:
1648 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1649 if res==None:
1650 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1651 else:
1652 if self.__numEquations==None: self.__numEquations=res[0]
1653 if self.__numSolutions==None: self.__numSolutions=res[1]
1654 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1655 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1656 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1657 for i,d in coefficients.iteritems():
1658 try:
1659 self.COEFFICIENTS[i].setValue(self.getDomain(),
1660 self.getNumEquations(),self.getNumSolutions(),
1661 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1662 self.alteredCoefficient(i)
1663 except IllegalCoefficientFunctionSpace,m:
1664 # if the function space is wrong then we try the reduced version:
1665 i_red=i+"_reduced"
1666 if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1667 try:
1668 self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1669 self.getNumEquations(),self.getNumSolutions(),
1670 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1671 self.alteredCoefficient(i_red)
1672 except IllegalCoefficientValue,m:
1673 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1674 except IllegalCoefficientFunctionSpace,m:
1675 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1676 else:
1677 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1678 except IllegalCoefficientValue,m:
1679 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1680 self.__altered_coefficients=True
1681 # check if the systrem is inhomogeneous:
1682 if len(coefficients)>0 and not self.isUsingLumping():
1683 q=self.getCoefficientOfGeneralPDE("q")
1684 r=self.getCoefficientOfGeneralPDE("r")
1685 homogeneous_constraint=True
1686 if not q.isEmpty() and not r.isEmpty():
1687 if util.Lsup(q*r)>0.:
1688 self.trace("Inhomogeneous constraint detected.")
1689 self.__invalidateSystem()
1690
1691 def getSystem(self):
1692 """
1693 return the operator and right hand side of the PDE
1694
1695 @return: the discrete version of the PDE
1696 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1697 """
1698 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1699 if self.isUsingLumping():
1700 if not self.__operator_is_Valid:
1701 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1702 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1703 if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1704 raise ValueError,"coefficient A in lumped matrix may not be present."
1705 if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1706 raise ValueError,"coefficient B in lumped matrix may not be present."
1707 if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1708 raise ValueError,"coefficient C in lumped matrix may not be present."
1709 if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1710 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1711 if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1712 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1713 if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1714 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1715 if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1716 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1717 if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1718 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1719 D=self.getCoefficientOfGeneralPDE("D")
1720 d=self.getCoefficientOfGeneralPDE("d")
1721 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1722 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1723 if not D.isEmpty():
1724 if self.getNumSolutions()>1:
1725 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1726 else:
1727 D_times_e=D
1728 else:
1729 D_times_e=escript.Data()
1730 if not d.isEmpty():
1731 if self.getNumSolutions()>1:
1732 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1733 else:
1734 d_times_e=d
1735 else:
1736 d_times_e=escript.Data()
1737
1738 if not D_reduced.isEmpty():
1739 if self.getNumSolutions()>1:
1740 D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1741 else:
1742 D_reduced_times_e=D_reduced
1743 else:
1744 D_reduced_times_e=escript.Data()
1745 if not d_reduced.isEmpty():
1746 if self.getNumSolutions()>1:
1747 d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1748 else:
1749 d_reduced_times_e=d_reduced
1750 else:
1751 d_reduced_times_e=escript.Data()
1752
1753 self.__operator=self.__getNewRightHandSide()
1754 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1755 self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1756 self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1757 else:
1758 self.getDomain().addPDEToRHS(self.__operator, \
1759 escript.Data(), \
1760 D_times_e, \
1761 d_times_e,\
1762 escript.Data())
1763 self.getDomain().addPDEToRHS(self.__operator, \
1764 escript.Data(), \
1765 D_reduced_times_e, \
1766 d_reduced_times_e,\
1767 escript.Data())
1768 self.__operator=1./self.__operator
1769 self.trace("New lumped operator has been built.")
1770 self.__operator_is_Valid=True
1771 if not self.__righthandside_isValid:
1772 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1773 self.getCoefficientOfGeneralPDE("X"), \
1774 self.getCoefficientOfGeneralPDE("Y"),\
1775 self.getCoefficientOfGeneralPDE("y"),\
1776 self.getCoefficientOfGeneralPDE("y_contact"))
1777 self.getDomain().addPDEToRHS(self.__righthandside, \
1778 self.getCoefficientOfGeneralPDE("X_reduced"), \
1779 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1780 self.getCoefficientOfGeneralPDE("y_reduced"),\
1781 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1782 self.trace("New right hand side as been built.")
1783 self.__righthandside_isValid=True
1784 else:
1785 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1786 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1787 self.getCoefficientOfGeneralPDE("A"), \
1788 self.getCoefficientOfGeneralPDE("B"), \
1789 self.getCoefficientOfGeneralPDE("C"), \
1790 self.getCoefficientOfGeneralPDE("D"), \
1791 self.getCoefficientOfGeneralPDE("X"), \
1792 self.getCoefficientOfGeneralPDE("Y"), \
1793 self.getCoefficientOfGeneralPDE("d"), \
1794 self.getCoefficientOfGeneralPDE("y"), \
1795 self.getCoefficientOfGeneralPDE("d_contact"), \
1796 self.getCoefficientOfGeneralPDE("y_contact"))
1797 self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
1798 self.getCoefficientOfGeneralPDE("A_reduced"), \
1799 self.getCoefficientOfGeneralPDE("B_reduced"), \
1800 self.getCoefficientOfGeneralPDE("C_reduced"), \
1801 self.getCoefficientOfGeneralPDE("D_reduced"), \
1802 self.getCoefficientOfGeneralPDE("X_reduced"), \
1803 self.getCoefficientOfGeneralPDE("Y_reduced"), \
1804 self.getCoefficientOfGeneralPDE("d_reduced"), \
1805 self.getCoefficientOfGeneralPDE("y_reduced"), \
1806 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1807 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1808 self.__applyConstraint()
1809 self.__righthandside=self.copyConstraint(self.__righthandside)
1810 self.trace("New system has been built.")
1811 self.__operator_is_Valid=True
1812 self.__righthandside_isValid=True
1813 elif not self.__righthandside_isValid:
1814 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1815 self.getCoefficientOfGeneralPDE("X"), \
1816 self.getCoefficientOfGeneralPDE("Y"),\
1817 self.getCoefficientOfGeneralPDE("y"),\
1818 self.getCoefficientOfGeneralPDE("y_contact"))
1819 self.getDomain().addPDEToRHS(self.__righthandside, \
1820 self.getCoefficientOfGeneralPDE("X_reduced"), \
1821 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1822 self.getCoefficientOfGeneralPDE("y_reduced"),\
1823 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1824 self.__righthandside=self.copyConstraint(self.__righthandside)
1825 self.trace("New right hand side has been built.")
1826 self.__righthandside_isValid=True
1827 elif not self.__operator_is_Valid:
1828 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1829 self.getCoefficientOfGeneralPDE("A"), \
1830 self.getCoefficientOfGeneralPDE("B"), \
1831 self.getCoefficientOfGeneralPDE("C"), \
1832 self.getCoefficientOfGeneralPDE("D"), \
1833 escript.Data(), \
1834 escript.Data(), \
1835 self.getCoefficientOfGeneralPDE("d"), \
1836 escript.Data(),\
1837 self.getCoefficientOfGeneralPDE("d_contact"), \
1838 escript.Data())
1839 self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \
1840 self.getCoefficientOfGeneralPDE("A_reduced"), \
1841 self.getCoefficientOfGeneralPDE("B_reduced"), \
1842 self.getCoefficientOfGeneralPDE("C_reduced"), \
1843 self.getCoefficientOfGeneralPDE("D_reduced"), \
1844 escript.Data(), \
1845 escript.Data(), \
1846 self.getCoefficientOfGeneralPDE("d_reduced"), \
1847 escript.Data(),\
1848 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1849 escript.Data())
1850 self.__applyConstraint()
1851 self.trace("New operator has been built.")
1852 self.__operator_is_Valid=True
1853 return (self.__operator,self.__righthandside)
1854
1855
1856 class Poisson(LinearPDE):
1857 """
1858 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1859
1860 M{-grad(grad(u)[j])[j] = f}
1861
1862 with natural boundary conditons
1863
1864 M{n[j]*grad(u)[j] = 0 }
1865
1866 and constraints:
1867
1868 M{u=0} where M{q>0}
1869
1870 """
1871
1872 def __init__(self,domain,debug=False):
1873 """
1874 initializes a new Poisson equation
1875
1876 @param domain: domain of the PDE
1877 @type domain: L{Domain<escript.Domain>}
1878 @param debug: if True debug informations are printed.
1879
1880 """
1881 super(Poisson, self).__init__(domain,1,1,debug)
1882 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1883 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1884 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1885 self.setSymmetryOn()
1886
1887 def setValue(self,**coefficients):
1888 """
1889 sets new values to coefficients
1890
1891 @param coefficients: new values assigned to coefficients
1892 @keyword f: value for right hand side M{f}
1893 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1894 @keyword q: mask for location of constraints
1895 @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1896 depending of reduced order is used for the representation of the equation.
1897 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1898 """
1899 super(Poisson, self).setValue(**coefficients)
1900
1901 def getCoefficientOfGeneralPDE(self,name):
1902 """
1903 return the value of the coefficient name of the general PDE
1904 @param name: name of the coefficient requested.
1905 @type name: C{string}
1906 @return: the value of the coefficient name
1907 @rtype: L{Data<escript.Data>}
1908 @raise IllegalCoefficient: if name is not one of coefficients
1909 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1910 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1911 """
1912 if name == "A" :
1913 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1914 elif name == "B" :
1915 return escript.Data()
1916 elif name == "C" :
1917 return escript.Data()
1918 elif name == "D" :
1919 return escript.Data()
1920 elif name == "X" :
1921 return escript.Data()
1922 elif name == "Y" :
1923 return self.getCoefficient("f")
1924 elif name == "d" :
1925 return escript.Data()
1926 elif name == "y" :
1927 return escript.Data()
1928 elif name == "d_contact" :
1929 return escript.Data()
1930 elif name == "y_contact" :
1931 return escript.Data()
1932 elif name == "A_reduced" :
1933 return escript.Data()
1934 elif name == "B_reduced" :
1935 return escript.Data()
1936 elif name == "C_reduced" :
1937 return escript.Data()
1938 elif name == "D_reduced" :
1939 return escript.Data()
1940 elif name == "X_reduced" :
1941 return escript.Data()
1942 elif name == "Y_reduced" :
1943 return self.getCoefficient("f_reduced")
1944 elif name == "d_reduced" :
1945 return escript.Data()
1946 elif name == "y_reduced" :
1947 return escript.Data()
1948 elif name == "d_contact_reduced" :
1949 return escript.Data()
1950 elif name == "y_contact_reduced" :
1951 return escript.Data()
1952 elif name == "r" :
1953 return escript.Data()
1954 elif name == "q" :
1955 return self.getCoefficient("q")
1956 else:
1957 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1958
1959 class Helmholtz(LinearPDE):
1960 """
1961 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1962
1963 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1964
1965 with natural boundary conditons
1966
1967 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1968
1969 and constraints:
1970
1971 M{u=r} where M{q>0}
1972
1973 """
1974
1975 def __init__(self,domain,debug=False):
1976 """
1977 initializes a new Poisson equation
1978
1979 @param domain: domain of the PDE
1980 @type domain: L{Domain<escript.Domain>}
1981 @param debug: if True debug informations are printed.
1982
1983 """
1984 super(Helmholtz, self).__init__(domain,1,1,debug)
1985 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1986 "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1987 "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1988 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1989 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1990 "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1991 "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1992 "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1993 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1994 self.setSymmetryOn()
1995
1996 def setValue(self,**coefficients):
1997 """
1998 sets new values to coefficients
1999
2000 @param coefficients: new values assigned to coefficients
2001 @keyword omega: value for coefficient M{S{omega}}
2002 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2003 @keyword k: value for coefficeint M{k}
2004 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2005 @keyword f: value for right hand side M{f}
2006 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2007 @keyword alpha: value for right hand side M{S{alpha}}
2008 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2009 @keyword g: value for right hand side M{g}
2010 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2011 @keyword r: prescribed values M{r} for the solution in constraints.
2012 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2013 depending of reduced order is used for the representation of the equation.
2014 @keyword q: mask for location of constraints
2015 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2016 depending of reduced order is used for the representation of the equation.
2017 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2018 """
2019 super(Helmholtz, self).setValue(**coefficients)
2020
2021 def getCoefficientOfGeneralPDE(self,name):
2022 """
2023 return the value of the coefficient name of the general PDE
2024
2025 @param name: name of the coefficient requested.
2026 @type name: C{string}
2027 @return: the value of the coefficient name
2028 @rtype: L{Data<escript.Data>}
2029 @raise IllegalCoefficient: if name is not one of coefficients
2030 "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2031 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2032 """
2033 if name == "A" :
2034 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2035 elif name == "B" :
2036 return escript.Data()
2037 elif name == "C" :
2038 return escript.Data()
2039 elif name == "D" :
2040 return self.getCoefficient("omega")
2041 elif name == "X" :
2042 return escript.Data()
2043 elif name == "Y" :
2044 return self.getCoefficient("f")
2045 elif name == "d" :
2046 return self.getCoefficient("alpha")
2047 elif name == "y" :
2048 return self.getCoefficient("g")
2049 elif name == "d_contact" :
2050 return escript.Data()
2051 elif name == "y_contact" :
2052 return escript.Data()
2053 elif name == "A_reduced" :
2054 return escript.Data()
2055 elif name == "B_reduced" :
2056 return escript.Data()
2057 elif name == "C_reduced" :
2058 return escript.Data()
2059 elif name == "D_reduced" :
2060 return escript.Data()
2061 elif name == "X_reduced" :
2062 return escript.Data()
2063 elif name == "Y_reduced" :
2064 return self.getCoefficient("f_reduced")
2065 elif name == "d_reduced" :
2066 return escript.Data()
2067 elif name == "y_reduced" :
2068 return self.getCoefficient("g_reduced")
2069 elif name == "d_contact_reduced" :
2070 return escript.Data()
2071 elif name == "y_contact_reduced" :
2072 return escript.Data()
2073 elif name == "r" :
2074 return self.getCoefficient("r")
2075 elif name == "q" :
2076 return self.getCoefficient("q")
2077 else:
2078 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2079
2080 class LameEquation(LinearPDE):
2081 """
2082 Class to define a Lame equation problem:
2083
2084 M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2085
2086 with natural boundary conditons:
2087
2088 M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
2089
2090 and constraints:
2091
2092 M{u[i]=r[i]} where M{q[i]>0}
2093
2094 """
2095
2096 def __init__(self,domain,debug=False):
2097 super(LameEquation, self).__init__(domain,\
2098 domain.getDim(),domain.getDim(),debug)
2099 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2100 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2101 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2102 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2103 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2104 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2105 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2106 self.setSymmetryOn()
2107
2108 def setValues(self,**coefficients):
2109 """
2110 sets new values to coefficients
2111
2112 @param coefficients: new values assigned to coefficients
2113 @keyword lame_mu: value for coefficient M{S{mu}}
2114 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2115 @keyword lame_lambda: value for coefficient M{S{lambda}}
2116 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2117 @keyword F: value for internal force M{F}
2118 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2119 @keyword sigma: value for initial stress M{S{sigma}}
2120 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2121 @keyword f: value for extrenal force M{f}
2122 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2123 @keyword r: prescribed values M{r} for the solution in constraints.
2124 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2125 depending of reduced order is used for the representation of the equation.
2126 @keyword q: mask for location of constraints
2127 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2128 depending of reduced order is used for the representation of the equation.
2129 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2130 """
2131 super(LameEquation, self).setValues(**coefficients)
2132
2133 def getCoefficientOfGeneralPDE(self,name):
2134 """
2135 return the value of the coefficient name of the general PDE
2136
2137 @param name: name of the coefficient requested.
2138 @type name: C{string}
2139 @return: the value of the coefficient name
2140 @rtype: L{Data<escript.Data>}
2141 @raise IllegalCoefficient: if name is not one of coefficients
2142 "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2143 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2144 """
2145 if name == "A" :
2146 out =self.createCoefficientOfGeneralPDE("A")
2147 for i in range(self.getDim()):
2148 for j in range(self.getDim()):
2149 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2150 out[i,j,j,i] += self.getCoefficient("lame_mu")
2151 out[i,j,i,j] += self.getCoefficient("lame_mu")
2152 return out
2153 elif name == "B" :
2154 return escript.Data()
2155 elif name == "C" :
2156 return escript.Data()
2157 elif name == "D" :
2158 return escript.Data()
2159 elif name == "X" :
2160 return self.getCoefficient("sigma")
2161 elif name == "Y" :
2162 return self.getCoefficient("F")
2163 elif name == "d" :
2164 return escript.Data()
2165 elif name == "y" :
2166 return self.getCoefficient("f")
2167 elif name == "d_contact" :
2168 return escript.Data()
2169 elif name == "y_contact" :
2170 return escript.Data()
2171 elif name == "A_reduced" :
2172 return escript.Data()
2173 elif name == "B_reduced" :
2174 return escript.Data()
2175 elif name == "C_reduced" :
2176 return escript.Data()
2177 elif name == "D_reduced" :
2178 return escript.Data()
2179 elif name == "X_reduced" :
2180 return escript.Data()
2181 elif name == "Y_reduced" :
2182 return escript.Data()
2183 elif name == "d_reduced" :
2184 return escript.Data()
2185 elif name == "y_reduced" :
2186 return escript.Data()
2187 elif name == "d_contact_reduced" :
2188 return escript.Data()
2189 elif name == "y_contact_reduced" :
2190 return escript.Data()
2191 elif name == "r" :
2192 return self.getCoefficient("r")
2193 elif name == "q" :
2194 return self.getCoefficient("q")
2195 else:
2196 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2197
2198 def LinearSinglePDE(domain,debug=False):
2199 """
2200 defines a single linear PDEs
2201
2202 @param domain: domain of the PDE
2203 @type domain: L{Domain<escript.Domain>}
2204 @param debug: if True debug informations are printed.
2205 @rtype: L{LinearPDE}
2206 """
2207 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2208
2209 def LinearPDESystem(domain,debug=False):
2210 """
2211 defines a system of linear PDEs
2212
2213 @param domain: domain of the PDE
2214 @type domain: L{Domain<escript.Domain>}
2215 @param debug: if True debug informations are printed.
2216 @rtype: L{LinearPDE}
2217 """
2218 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2219
2220 class TransportPDE(object):
2221 """
2222 Warning: This is still a very experimental. The class is still changing!
2223
2224 Mu_{,t} =-(A_{ij}u_{,j})_j-(B_{j}u)_{,j} + C_{j} u_{,j} + Y_i + X_{i,i}
2225
2226 u=r where q>0
2227
2228 all coefficients are constant over time.
2229
2230 typical usage:
2231
2232 p=TransportPDE(dom)
2233 p.setValue(M=Scalar(1.,Function(dom),C=Scalar(1.,Function(dom)*[-1.,0.])
2234 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2235 t=0
2236 dt=0.1
2237 while (t<1.):
2238 u=p.solve(dt)
2239
2240 """
2241 def __init__(self,domain,num_equations=1,theta=0.5,useSUPG=False,trace=True):
2242 self.__domain=domain
2243 self.__num_equations=num_equations
2244 self.__useSUPG=useSUPG
2245 self.__trace=trace
2246 self.__theta=theta
2247 self.__matrix_type=0
2248 self.__reduced=True
2249 self.__reassemble=True
2250 if self.__useSUPG:
2251 self.__pde=LinearPDE(domain,numEquations=num_equations,numSolutions=num_equations,debug=trace)
2252 self.__pde.setSymmetryOn()
2253 self.__pde.setReducedOrderOn()
2254 else:
2255 self.__transport_problem=self.__getNewTransportProblem()
2256 self.setTolerance()
2257 self.__M=escript.Data()
2258 self.__A=escript.Data()
2259 self.__B=escript.Data()
2260 self.__C=escript.Data()
2261 self.__D=escript.Data()
2262 self.__X=escript.Data()
2263 self.__Y=escript.Data()
2264 self.__d=escript.Data()
2265 self.__y=escript.Data()
2266 self.__d_contact=escript.Data()
2267 self.__y_contact=escript.Data()
2268 self.__r=escript.Data()
2269 self.__q=escript.Data()
2270
2271 def trace(self,text):
2272 if self.__trace: print text
2273 def getSafeTimeStepSize(self):
2274 if self.__useSUPG:
2275 if self.__reassemble:
2276 h=self.__domain.getSize()
2277 dt=None
2278 if not self.__A.isEmpty():
2279 dt2=util.inf(h**2*self.__M/util.length(self.__A))
2280 if dt == None:
2281 dt = dt2
2282 else:
2283 dt=1./(1./dt+1./dt2)
2284 if not self.__B.isEmpty():
2285 dt2=util.inf(h*self.__M/util.length(self.__B))
2286 if dt == None:
2287 dt = dt2
2288 else:
2289 dt=1./(1./dt+1./dt2)
2290 if not self.__C.isEmpty():
2291 dt2=util.inf(h*self.__M/util.length(self.__C))
2292 if dt == None:
2293 dt = dt2
2294 else:
2295 dt=1./(1./dt+1./dt2)
2296 if not self.__D.isEmpty():
2297 dt2=util.inf(self.__M/util.length(self.__D))
2298 if dt == None:
2299 dt = dt2
2300 else:
2301 dt=1./(1./dt+1./dt2)
2302 self.__dt = dt/2
2303 return self.__dt
2304 else:
2305 return self.__getTransportProblem().getSafeTimeStepSize()
2306 def getDomain(self):
2307 return self.__domain
2308 def getTheta(self):
2309 return self.__theta
2310 def getNumEquations(self):
2311 return self.__num_equations
2312 def setReducedOn(self):
2313 if not self.reduced():
2314 if self.__useSUPG:
2315 self.__pde.setReducedOrderOn()
2316 else:
2317 self.__transport_problem=self.__getNewTransportProblem()
2318 self.__reduced=True
2319 def setReducedOff(self):
2320 if self.reduced():
2321 if self.__useSUPG:
2322 self.__pde.setReducedOrderOff()
2323 else:
2324 self.__transport_problem=self.__getNewTransportProblem()
2325 self.__reduced=False
2326 def reduced(self):
2327 return self.__reduced
2328 def getFunctionSpace(self):
2329 if self.reduced():
2330 return escript.ReducedSolution(self.getDomain())
2331 else:
2332 return escript.Solution(self.getDomain())
2333
2334 def setTolerance(self,tol=1.e-8):
2335 self.__tolerance=tol
2336 if self.__useSUPG:
2337 self.__pde.setTolerance(self.__tolerance)
2338
2339 def __getNewTransportProblem(self):
2340 """
2341 returns an instance of a new operator
2342 """
2343 self.trace("New Transport problem is allocated.")
2344 return self.getDomain().newTransportProblem( \
2345 self.getTheta(),
2346 self.getNumEquations(), \
2347 self.getFunctionSpace(), \
2348 self.__matrix_type)
2349
2350 def __getNewSolutionVector(self):
2351 if self.getNumEquations() ==1 :
2352 out=escript.Data(0.0,(),self.getFunctionSpace())
2353 else:
2354 out=escript.Data(0.0,(self.getNumEquations(),),self.getFunctionSpace())
2355 return out
2356
2357 def __getTransportProblem(self):
2358 if self.__reassemble:
2359 self.__source=self.__getNewSolutionVector()
2360 self.__transport_problem.reset()
2361 self.getDomain().addPDEToTransportProblem(
2362 self.__transport_problem,
2363 self.__source,
2364 self.__M,
2365 self.__A,
2366 self.__B,
2367 self.__C,
2368 self.__D,
2369 self.__X,
2370 self.__Y,
2371 self.__d,
2372 self.__y,
2373 self.__d_contact,
2374 self.__y_contact)
2375 self.__transport_problem.insertConstraint(self.__source,self.__q,self.__r)
2376 self.__reassemble=False
2377 return self.__transport_problem
2378 def setValue(self,M=None, A=None, B=None, C=None, D=None, X=None, Y=None,
2379 d=None, y=None, d_contact=None, y_contact=None, q=None, r=None):
2380 if not M==None:
2381 self.__reassemble=True
2382 self.__M=M
2383 if not A==None:
2384 self.__reassemble=True
2385 self.__A=A
2386 if not B==None:
2387 self.__reassemble=True
2388 self.__B=B
2389 if not C==None:
2390 self.__reassemble=True
2391 self.__C=C
2392 if not D==None:
2393 self.__reassemble=True
2394 self.__D=D
2395 if not X==None:
2396 self.__reassemble=True
2397 self.__X=X
2398 if not Y==None:
2399 self.__reassemble=True
2400 self.__Y=Y
2401 if not d==None:
2402 self.__reassemble=True
2403 self.__d=d
2404 if not y==None:
2405 self.__reassemble=True
2406 self.__y=y
2407 if not d_contact==None:
2408 self.__reassemble=True
2409 self.__d_contact=d_contact
2410 if not y_contact==None:
2411 self.__reassemble=True
2412 self.__y_contact=y_contact
2413 if not q==None:
2414 self.__reassemble=True
2415 self.__q=q
2416 if not r==None:
2417 self.__reassemble=True
2418 self.__r=r
2419
2420 def setInitialSolution(self,u):
2421 if self.__useSUPG:
2422 self.__u=util.interpolate(u,self.getFunctionSpace())
2423 else:
2424 self.__transport_problem.setInitialValue(util.interpolate(u,self.getFunctionSpace()))
2425
2426 def solve(self,dt,**kwarg):
2427 if self.__useSUPG:
2428 if self.__reassemble:
2429 self.__pde.setValue(D=self.__M,d=self.__d,d_contact=self.__d_contact,q=self.__q) # ,r=self.__r)
2430 self.__reassemble=False
2431 dt2=self.getSafeTimeStepSize()
2432 nn=max(math.ceil(dt/self.getSafeTimeStepSize()),1.)
2433 dt2=dt/nn
2434 nnn=0
2435 u=self.__u
2436 self.trace("number of substeps is %d."%nn)
2437 while nnn<nn :
2438 self.__setSUPG(u,u,dt2/2)
2439 u_half=self.__pde.getSolution(verbose=True)
2440 self.__setSUPG(u,u_half,dt2)
2441 u=self.__pde.getSolution(verbose=True)
2442 nnn+=1
2443 self.__u=u
2444 return self.__u
2445 else:
2446 kwarg["tolerance"]=self.__tolerance
2447 tp=self.__getTransportProblem()
2448 return tp.solve(self.__source,dt,kwarg)
2449 def __setSUPG(self,u0,u,dt):
2450 g=util.grad(u)
2451 X=0
2452 Y=self.__M*u0
2453 X=0
2454 self.__pde.setValue(r=u0)
2455 if not self.__A.isEmpty():
2456 X=X+dt*util.matrixmult(self.__A,g)
2457 if not self.__B.isEmpty():
2458 X=X+dt*self.__B*u
2459 if not self.__C.isEmpty():
2460 Y=Y+dt*util.inner(self.__C,g)
2461 if not self.__D.isEmpty():
2462 Y=Y+dt*self.__D*u
2463 if not self.__X.isEmpty():
2464 X=X+dt*self.__X
2465 if not self.__Y.isEmpty():
2466 Y=Y+dt*self.__Y
2467 self.__pde.setValue(X=X,Y=Y)
2468 if not self.__y.isEmpty():
2469 self.__pde.setValue(y=dt*self.__y)
2470 if not self.__y_contact.isEmpty():
2471 self.__pde.setValue(y=dt*self.__y_contact)
2472 self.__pde.setValue(r=u0)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26