/[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 1639 - (show annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 112746 byte(s)


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