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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26