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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26