/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1819 - (show annotations)
Tue Sep 30 05:58:06 2008 UTC (11 years, 1 month ago) by artak
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 113062 byte(s)
Firs version of symmetric Gauss-Seidel preconditioner with coloring
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 """
23 The module provides an interface to define and solve linear partial
24 differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
25 solver capabilities in itself but hands the PDE over to
26 the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
27 The general interface is provided through the L{LinearPDE} class. The
28 L{AdvectivePDE} which is derived from the L{LinearPDE} class
29 provides an interface to PDE dominated by its advective terms. The L{Poisson},
30 L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
31 classs which are also derived form the L{LinearPDE} class should be used
32 to define of solve these sepecial PDEs.
33
34 @var __author__: name of author
35 @var __copyright__: copyrights
36 @var __license__: licence agreement
37 @var __url__: url entry point on documentation
38 @var __version__: version
39 @var __date__: date of the version
40 """
41
42 import math
43 import escript
44 import util
45 import numarray
46
47 __author__="Lutz Gross, l.gross@uq.edu.au"
48
49
50 class IllegalCoefficient(ValueError):
51 """
52 raised if an illegal coefficient of the general ar particular PDE is requested.
53 """
54 pass
55
56 class IllegalCoefficientValue(ValueError):
57 """
58 raised if an incorrect value for a coefficient is used.
59 """
60 pass
61
62 class IllegalCoefficientFunctionSpace(ValueError):
63 """
64 raised if an incorrect function space for a coefficient is used.
65 """
66
67 class UndefinedPDEError(ValueError):
68 """
69 raised if a PDE is not fully defined yet.
70 """
71 pass
72
73 class PDECoefficient(object):
74 """
75 A class for describing a PDE coefficient
76
77 @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
78 @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
79 @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
80 @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
81 @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
82 @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
83 @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
84 @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
85 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
86 @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
87 @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
88 @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
89 @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
90 @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
91
92 """
93 INTERIOR=0
94 BOUNDARY=1
95 CONTACT=2
96 SOLUTION=3
97 REDUCED=4
98 BY_EQUATION=5
99 BY_SOLUTION=6
100 BY_DIM=7
101 OPERATOR=10
102 RIGHTHANDSIDE=11
103 BOTH=12
104 INTERIOR_REDUCED=13
105 BOUNDARY_REDUCED=14
106 CONTACT_REDUCED=15
107
108 def __init__(self, where, pattern, altering):
109 """
110 Initialise a PDE Coefficient type
111
112 @param where: describes where the coefficient lives
113 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
114 L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
115 @param pattern: describes the shape of the coefficient and how the shape is build for a given
116 spatial dimension and numbers of equation and solution in then PDE. For instance,
117 (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
118 is instanciated as shape (3,2,2) in case of a three equations and two solution components
119 on a 2-dimensional domain. In the case of single equation and a single solution component
120 the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
121 the example would be read as (2,).
122 @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
123 @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
124 @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
125 @param reduced: indicates if reduced
126 @type reduced: C{bool}
127 """
128 super(PDECoefficient, self).__init__()
129 self.what=where
130 self.pattern=pattern
131 self.altering=altering
132 self.resetValue()
133
134 def resetValue(self):
135 """
136 resets coefficient value to default
137 """
138 self.value=escript.Data()
139
140 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
141 """
142 defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
143
144 @param domain: domain on which the PDE uses the coefficient
145 @type domain: L{Domain<escript.Domain>}
146 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
147 @type reducedEquationOrder: C{bool}
148 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
149 @type reducedSolutionOrder: C{bool}
150 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
151 @rtype: L{FunctionSpace<escript.FunctionSpace>}
152 """
153 if self.what==self.INTERIOR:
154 return escript.Function(domain)
155 elif self.what==self.INTERIOR_REDUCED:
156 return escript.ReducedFunction(domain)
157 elif self.what==self.BOUNDARY:
158 return escript.FunctionOnBoundary(domain)
159 elif self.what==self.BOUNDARY_REDUCED:
160 return escript.ReducedFunctionOnBoundary(domain)
161 elif self.what==self.CONTACT:
162 return escript.FunctionOnContactZero(domain)
163 elif self.what==self.CONTACT_REDUCED:
164 return escript.ReducedFunctionOnContactZero(domain)
165 elif self.what==self.SOLUTION:
166 if reducedEquationOrder and reducedSolutionOrder:
167 return escript.ReducedSolution(domain)
168 else:
169 return escript.Solution(domain)
170 elif self.what==self.REDUCED:
171 return escript.ReducedSolution(domain)
172
173 def getValue(self):
174 """
175 returns the value of the coefficient
176
177 @return: value of the coefficient
178 @rtype: L{Data<escript.Data>}
179 """
180 return self.value
181
182 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
183 """
184 set the value of the coefficient to a new value
185
186 @param domain: domain on which the PDE uses the coefficient
187 @type domain: L{Domain<escript.Domain>}
188 @param numEquations: number of equations of the PDE
189 @type numEquations: C{int}
190 @param numSolutions: number of components of the PDE solution
191 @type numSolutions: C{int}
192 @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
193 @type reducedEquationOrder: C{bool}
194 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
195 @type reducedSolutionOrder: C{bool}
196 @param newValue: number of components of the PDE solution
197 @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
198 @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
199 @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
200 """
201 if newValue==None:
202 newValue=escript.Data()
203 elif isinstance(newValue,escript.Data):
204 if not newValue.isEmpty():
205 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
206 try:
207 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
208 except:
209 raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
210 else:
211 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
212 if not newValue.isEmpty():
213 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
214 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
215 self.value=newValue
216
217 def isAlteringOperator(self):
218 """
219 checks if the coefficient alters the operator of the PDE
220
221 @return: True if the operator of the PDE is changed when the coefficient is changed
222 @rtype: C{bool}
223 """
224 if self.altering==self.OPERATOR or self.altering==self.BOTH:
225 return not None
226 else:
227 return None
228
229 def isAlteringRightHandSide(self):
230 """
231 checks if the coefficeint alters the right hand side of the PDE
232
233 @rtype: C{bool}
234 @return: True if the right hand side of the PDE is changed when the coefficient is changed
235 """
236 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
237 return not None
238 else:
239 return None
240
241 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
242 """
243 tries to estimate the number of equations and number of solutions if the coefficient has the given shape
244
245 @param domain: domain on which the PDE uses the coefficient
246 @type domain: L{Domain<escript.Domain>}
247 @param shape: suggested shape of the coefficient
248 @type shape: C{tuple} of C{int} values
249 @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
250 If no appropriate numbers could be identified, C{None} is returned
251 @rtype: C{tuple} of two C{int} values or C{None}
252 """
253 dim=domain.getDim()
254 if len(shape)>0:
255 num=max(shape)+1
256 else:
257 num=1
258 search=[]
259 if self.definesNumEquation() and self.definesNumSolutions():
260 for u in range(num):
261 for e in range(num):
262 search.append((e,u))
263 search.sort(self.__CompTuple2)
264 for item in search:
265 s=self.getShape(domain,item[0],item[1])
266 if len(s)==0 and len(shape)==0:
267 return (1,1)
268 else:
269 if s==shape: return item
270 elif self.definesNumEquation():
271 for e in range(num,0,-1):
272 s=self.getShape(domain,e,0)
273 if len(s)==0 and len(shape)==0:
274 return (1,None)
275 else:
276 if s==shape: return (e,None)
277
278 elif self.definesNumSolutions():
279 for u in range(num,0,-1):
280 s=self.getShape(domain,0,u)
281 if len(s)==0 and len(shape)==0:
282 return (None,1)
283 else:
284 if s==shape: return (None,u)
285 return None
286 def definesNumSolutions(self):
287 """
288 checks if the coefficient allows to estimate the number of solution components
289
290 @return: True if the coefficient allows an estimate of the number of solution components
291 @rtype: C{bool}
292 """
293 for i in self.pattern:
294 if i==self.BY_SOLUTION: return True
295 return False
296
297 def definesNumEquation(self):
298 """
299 checks if the coefficient allows to estimate the number of equations
300
301 @return: True if the coefficient allows an estimate of the number of equations
302 @rtype: C{bool}
303 """
304 for i in self.pattern:
305 if i==self.BY_EQUATION: return True
306 return False
307
308 def __CompTuple2(self,t1,t2):
309 """
310 Compare two tuples of possible number of equations and number of solutions
311
312 @param t1: The first tuple
313 @param t2: The second tuple
314
315 """
316
317 dif=t1[0]+t1[1]-(t2[0]+t2[1])
318 if dif<0: return 1
319 elif dif>0: return -1
320 else: return 0
321
322 def getShape(self,domain,numEquations=1,numSolutions=1):
323 """
324 builds the required shape of the coefficient
325
326 @param domain: domain on which the PDE uses the coefficient
327 @type domain: L{Domain<escript.Domain>}
328 @param numEquations: number of equations of the PDE
329 @type numEquations: C{int}
330 @param numSolutions: number of components of the PDE solution
331 @type numSolutions: C{int}
332 @return: shape of the coefficient
333 @rtype: C{tuple} of C{int} values
334 """
335 dim=domain.getDim()
336 s=()
337 for i in self.pattern:
338 if i==self.BY_EQUATION:
339 if numEquations>1: s=s+(numEquations,)
340 elif i==self.BY_SOLUTION:
341 if numSolutions>1: s=s+(numSolutions,)
342 else:
343 s=s+(dim,)
344 return s
345
346 class LinearPDE(object):
347 """
348 This class is used to define a general linear, steady, second order PDE
349 for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
350
351 For a single PDE with a solution with a single component the linear PDE is defined in the following form:
352
353 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)}
354
355
356 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
357 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
358 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
359 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
360 L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
361 such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B_reduced}, M{C_reduced}, M{X_reduced}
362 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.
363
364 The following natural boundary conditions are considered:
365
366 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}
367
368 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>}.
369
370
371 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
372
373 M{u=r} where M{q>0}
374
375 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
376 The constraints override any other condition set by the PDE or the boundary condition.
377
378 The PDE is symmetrical if
379
380 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]}
381
382 For a system of PDEs and a solution with several components the PDE has the form
383
384 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] }
385
386 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.
387 The natural boundary conditions take the form:
388
389 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]}
390
391
392 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>}.
393
394 Constraints take the form
395
396 M{u[i]=r[i]} where M{q[i]>0}
397
398 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
399
400 The system of PDEs is symmetrical if
401
402 - M{A[i,j,k,l]=A[k,l,i,j]}
403 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
404 - M{B[i,j,k]=C[k,i,j]}
405 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
406 - M{D[i,k]=D[i,k]}
407 - M{D_reduced[i,k]=D_reduced[i,k]}
408 - M{d[i,k]=d[k,i]}
409 - M{d_reduced[i,k]=d_reduced[k,i]}
410
411 L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
412 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
413 defined as
414
415 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]}
416
417 For the case of single solution component and single PDE M{J} is defined
418
419 M{J_{j}=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
420
421 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
422 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
423 the contact condition takes the form
424
425 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]}
426
427 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
428 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
429 L{jump<util.jump>}.
430 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>}.
431 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>}.
432 In case of a single PDE and a single component solution the contact condition takes the form
433
434 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
435
436 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>}
437
438 @cvar DEFAULT: The default method used to solve the system of linear equations
439 @cvar DIRECT: The direct solver based on LDU factorization
440 @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
441 @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
442 @cvar CR: The conjugate residual method
443 @cvar CGS: The conjugate gardient square method
444 @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
445 @cvar TFQMR: Transport Free Quasi Minimal Residual method.
446 @cvar MINRES: Minimum residual method.
447 @cvar SSOR: The symmetric overrealaxtion method
448 @cvar ILU0: The incomplete LU factorization preconditioner with no fill in
449 @cvar ILUT: The incomplete LU factorization preconditioner with will in
450 @cvar JACOBI: The Jacobi preconditioner
451 @cvar GMRES: The Gram-Schmidt minimum residual method
452 @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
453 @cvar LUMPING: Matrix lumping.
454 @cvar NO_REORDERING: No matrix reordering allowed
455 @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
456 @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
457 @cvar PASO: PASO solver package
458 @cvar SCSL: SGI SCSL solver library
459 @cvar MKL: Intel's MKL solver library
460 @cvar UMFPACK: the UMFPACK library
461 @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
462 @cvar ITERATIVE: The default iterative solver
463 @cvar AMG: algebraic multi grid
464 @cvar RILU: recursive ILU
465 @cvar GS: Gauss-Seidel solver
466
467 """
468 DEFAULT= 0
469 DIRECT= 1
470 CHOLEVSKY= 2
471 PCG= 3
472 CR= 4
473 CGS= 5
474 BICGSTAB= 6
475 SSOR= 7
476 ILU0= 8
477 ILUT= 9
478 JACOBI= 10
479 GMRES= 11
480 PRES20= 12
481 LUMPING= 13
482 NO_REORDERING= 17
483 MINIMUM_FILL_IN= 18
484 NESTED_DISSECTION= 19
485 SCSL= 14
486 MKL= 15
487 UMFPACK= 16
488 ITERATIVE= 20
489 PASO= 21
490 AMG= 22
491 RILU = 23
492 TRILINOS = 24
493 NONLINEAR_GMRES = 25
494 TFQMR = 26
495 MINRES = 27
496 GS=28
497
498 SMALL_TOLERANCE=1.e-13
499 __PACKAGE_KEY="package"
500 __METHOD_KEY="method"
501 __SYMMETRY_KEY="symmetric"
502 __TOLERANCE_KEY="tolerance"
503 __PRECONDITIONER_KEY="preconditioner"
504
505
506 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
507 """
508 initializes a new linear PDE
509
510 @param domain: domain of the PDE
511 @type domain: L{Domain<escript.Domain>}
512 @param numEquations: number of equations. If numEquations==None the number of equations
513 is exracted from the PDE coefficients.
514 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
515 is exracted from the PDE coefficients.
516 @param debug: if True debug informations are printed.
517
518 """
519 super(LinearPDE, self).__init__()
520 #
521 # the coefficients of the general PDE:
522 #
523 self.__COEFFICIENTS_OF_GENEARL_PDE={
524 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
525 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
526 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
527 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
528 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
529 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
530 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
531 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
532 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
533 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
534 "A_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
535 "B_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
536 "C_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
537 "D_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
538 "X_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
539 "Y_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
540 "d_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
541 "y_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
542 "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
543 "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
544 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
545 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
546
547 # COEFFICIENTS can be overwritten by subclasses:
548 self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
549 self.__altered_coefficients=False
550 # initialize attributes
551 self.__debug=debug
552 self.__domain=domain
553 self.__numEquations=numEquations
554 self.__numSolutions=numSolutions
555 self.__resetSystem()
556
557 # set some default values:
558 self.__reduce_equation_order=False
559 self.__reduce_solution_order=False
560 self.__tolerance=1.e-8
561 self.__solver_method=self.DEFAULT
562 self.__solver_package=self.DEFAULT
563 self.__preconditioner=self.DEFAULT
564 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
565 self.__sym=False
566
567 self.resetCoefficients()
568 self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
569 # =============================================================================
570 # general stuff:
571 # =============================================================================
572 def __str__(self):
573 """
574 returns string representation of the PDE
575
576 @return: a simple representation of the PDE
577 @rtype: C{str}
578 """
579 return "<LinearPDE %d>"%id(self)
580 # =============================================================================
581 # debug :
582 # =============================================================================
583 def setDebugOn(self):
584 """
585 switches on debugging
586 """
587 self.__debug=not None
588
589 def setDebugOff(self):
590 """
591 switches off debugging
592 """
593 self.__debug=None
594
595 def trace(self,text):
596 """
597 print the text message if debugging is swiched on.
598 @param text: message
599 @type text: C{string}
600 """
601 if self.__debug: print "%s: %s"%(str(self),text)
602
603 # =============================================================================
604 # some service functions:
605 # =============================================================================
606 def getDomain(self):
607 """
608 returns the domain of the PDE
609
610 @return: the domain of the PDE
611 @rtype: L{Domain<escript.Domain>}
612 """
613 return self.__domain
614
615 def getDim(self):
616 """
617 returns the spatial dimension of the PDE
618
619 @return: the spatial dimension of the PDE domain
620 @rtype: C{int}
621 """
622 return self.getDomain().getDim()
623
624 def getNumEquations(self):
625 """
626 returns the number of equations
627
628 @return: the number of equations
629 @rtype: C{int}
630 @raise UndefinedPDEError: if the number of equations is not be specified yet.
631 """
632 if self.__numEquations==None:
633 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
634 else:
635 return self.__numEquations
636
637 def getNumSolutions(self):
638 """
639 returns the number of unknowns
640
641 @return: the number of unknowns
642 @rtype: C{int}
643 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
644 """
645 if self.__numSolutions==None:
646 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
647 else:
648 return self.__numSolutions
649
650 def reduceEquationOrder(self):
651 """
652 return status for order reduction for equation
653
654 @return: return True is reduced interpolation order is used for the represenation of the equation
655 @rtype: L{bool}
656 """
657 return self.__reduce_equation_order
658
659 def reduceSolutionOrder(self):
660 """
661 return status for order reduction for the solution
662
663 @return: return True is reduced interpolation order is used for the represenation of the solution
664 @rtype: L{bool}
665 """
666 return self.__reduce_solution_order
667
668 def getFunctionSpaceForEquation(self):
669 """
670 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
671
672 @return: representation space of equation
673 @rtype: L{FunctionSpace<escript.FunctionSpace>}
674 """
675 if self.reduceEquationOrder():
676 return escript.ReducedSolution(self.getDomain())
677 else:
678 return escript.Solution(self.getDomain())
679
680 def getFunctionSpaceForSolution(self):
681 """
682 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
683
684 @return: representation space of solution
685 @rtype: L{FunctionSpace<escript.FunctionSpace>}
686 """
687 if self.reduceSolutionOrder():
688 return escript.ReducedSolution(self.getDomain())
689 else:
690 return escript.Solution(self.getDomain())
691
692
693 def getOperator(self):
694 """
695 provides access to the operator of the PDE
696
697 @return: the operator of the PDE
698 @rtype: L{Operator<escript.Operator>}
699 """
700 m=self.getSystem()[0]
701 if self.isUsingLumping():
702 return self.copyConstraint(1./m)
703 else:
704 return m
705
706 def getRightHandSide(self):
707 """
708 provides access to the right hand side of the PDE
709 @return: the right hand side of the PDE
710 @rtype: L{Data<escript.Data>}
711 """
712 r=self.getSystem()[1]
713 if self.isUsingLumping():
714 return self.copyConstraint(r)
715 else:
716 return r
717
718 def applyOperator(self,u=None):
719 """
720 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
721
722 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
723 the current solution is used.
724 @type u: L{Data<escript.Data>} or None
725 @return: image of u
726 @rtype: L{Data<escript.Data>}
727 """
728 if u==None:
729 return self.getOperator()*self.getSolution()
730 else:
731 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
732
733 def getResidual(self,u=None):
734 """
735 return the residual of u or the current solution if u is not present.
736
737 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
738 the current solution is used.
739 @type u: L{Data<escript.Data>} or None
740 @return: residual of u
741 @rtype: L{Data<escript.Data>}
742 """
743 return self.applyOperator(u)-self.getRightHandSide()
744
745 def checkSymmetry(self,verbose=True):
746 """
747 test the PDE for symmetry.
748
749 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
750 @type verbose: C{bool}
751 @return: True if the PDE is symmetric.
752 @rtype: L{Data<escript.Data>}
753 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
754 """
755 verbose=verbose or self.__debug
756 out=True
757 if self.getNumSolutions()!=self.getNumEquations():
758 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
759 out=False
760 else:
761 A=self.getCoefficientOfGeneralPDE("A")
762 if not A.isEmpty():
763 tol=util.Lsup(A)*self.SMALL_TOLERANCE
764 if self.getNumSolutions()>1:
765 for i in range(self.getNumEquations()):
766 for j in range(self.getDim()):
767 for k in range(self.getNumSolutions()):
768 for l in range(self.getDim()):
769 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
770 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)
771 out=False
772 else:
773 for j in range(self.getDim()):
774 for l in range(self.getDim()):
775 if util.Lsup(A[j,l]-A[l,j])>tol:
776 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
777 out=False
778 B=self.getCoefficientOfGeneralPDE("B")
779 C=self.getCoefficientOfGeneralPDE("C")
780 if B.isEmpty() and not C.isEmpty():
781 if verbose: print "non-symmetric PDE because B is not present but C is"
782 out=False
783 elif not B.isEmpty() and C.isEmpty():
784 if verbose: print "non-symmetric PDE because C is not present but B is"
785 out=False
786 elif not B.isEmpty() and not C.isEmpty():
787 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
788 if self.getNumSolutions()>1:
789 for i in range(self.getNumEquations()):
790 for j in range(self.getDim()):
791 for k in range(self.getNumSolutions()):
792 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
793 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
794 out=False
795 else:
796 for j in range(self.getDim()):
797 if util.Lsup(B[j]-C[j])>tol:
798 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
799 out=False
800 if self.getNumSolutions()>1:
801 D=self.getCoefficientOfGeneralPDE("D")
802 if not D.isEmpty():
803 tol=util.Lsup(D)*self.SMALL_TOLERANCE
804 for i in range(self.getNumEquations()):
805 for k in range(self.getNumSolutions()):
806 if util.Lsup(D[i,k]-D[k,i])>tol:
807 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
808 out=False
809 d=self.getCoefficientOfGeneralPDE("d")
810 if not d.isEmpty():
811 tol=util.Lsup(d)*self.SMALL_TOLERANCE
812 for i in range(self.getNumEquations()):
813 for k in range(self.getNumSolutions()):
814 if util.Lsup(d[i,k]-d[k,i])>tol:
815 if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
816 out=False
817 d_contact=self.getCoefficientOfGeneralPDE("d_contact")
818 if not d_contact.isEmpty():
819 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
820 for i in range(self.getNumEquations()):
821 for k in range(self.getNumSolutions()):
822 if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
823 if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
824 out=False
825 # and now the reduced coefficients
826 A_reduced=self.getCoefficientOfGeneralPDE("A_reduced")
827 if not A_reduced.isEmpty():
828 tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE
829 if self.getNumSolutions()>1:
830 for i in range(self.getNumEquations()):
831 for j in range(self.getDim()):
832 for k in range(self.getNumSolutions()):
833 for l in range(self.getDim()):
834 if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol:
835 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)
836 out=False
837 else:
838 for j in range(self.getDim()):
839 for l in range(self.getDim()):
840 if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol:
841 if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j)
842 out=False
843 B_reduced=self.getCoefficientOfGeneralPDE("B_reduced")
844 C_reduced=self.getCoefficientOfGeneralPDE("C_reduced")
845 if B_reduced.isEmpty() and not C_reduced.isEmpty():
846 if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is"
847 out=False
848 elif not B_reduced.isEmpty() and C_reduced.isEmpty():
849 if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is"
850 out=False
851 elif not B_reduced.isEmpty() and not C_reduced.isEmpty():
852 tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2.
853 if self.getNumSolutions()>1:
854 for i in range(self.getNumEquations()):
855 for j in range(self.getDim()):
856 for k in range(self.getNumSolutions()):
857 if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol:
858 if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j)
859 out=False
860 else:
861 for j in range(self.getDim()):
862 if util.Lsup(B_reduced[j]-C_reduced[j])>tol:
863 if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j)
864 out=False
865 if self.getNumSolutions()>1:
866 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
867 if not D_reduced.isEmpty():
868 tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE
869 for i in range(self.getNumEquations()):
870 for k in range(self.getNumSolutions()):
871 if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol:
872 if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i)
873 out=False
874 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
875 if not d_reduced.isEmpty():
876 tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE
877 for i in range(self.getNumEquations()):
878 for k in range(self.getNumSolutions()):
879 if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol:
880 if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i)
881 out=False
882 d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
883 if not d_contact_reduced.isEmpty():
884 tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE
885 for i in range(self.getNumEquations()):
886 for k in range(self.getNumSolutions()):
887 if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol:
888 if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i)
889 out=False
890 return out
891
892 def getSolution(self,**options):
893 """
894 returns the solution of the PDE. If the solution is not valid the PDE is solved.
895
896 @return: the solution
897 @rtype: L{Data<escript.Data>}
898 @param options: solver options
899 @keyword verbose: True to get some information during PDE solution
900 @type verbose: C{bool}
901 @keyword reordering: reordering scheme to be used during elimination. Allowed values are
902 L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
903 @keyword iter_max: maximum number of iteration steps allowed.
904 @keyword drop_tolerance: threshold for drupping in L{ILUT}
905 @keyword drop_storage: maximum of allowed memory in L{ILUT}
906 @keyword truncation: maximum number of residuals in L{GMRES}
907 @keyword restart: restart cycle length in L{GMRES}
908 """
909 if not self.__solution_isValid:
910 mat,f=self.getSystem()
911 if self.isUsingLumping():
912 self.__solution=self.copyConstraint(f*mat)
913 else:
914 options[self.__TOLERANCE_KEY]=self.getTolerance()
915 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
916 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
917 options[self.__PACKAGE_KEY]=self.getSolverPackage()
918 options[self.__SYMMETRY_KEY]=self.isSymmetric()
919 self.trace("PDE is resolved.")
920 self.trace("solver options: %s"%str(options))
921 self.__solution=mat.solve(f,options)
922 self.__solution_isValid=True
923 return self.__solution
924
925 def getFlux(self,u=None):
926 """
927 returns the flux M{J} for a given M{u}
928
929 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]}
930
931 or
932
933 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
934
935 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
936 @type u: L{Data<escript.Data>} or None
937 @return: flux
938 @rtype: L{Data<escript.Data>}
939 """
940 if u==None: u=self.getSolution()
941 return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \
942 +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \
943 -util.self.getCoefficientOfGeneralPDE("X") \
944 +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
945 +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \
946 -util.self.getCoefficientOfGeneralPDE("X_reduced")
947 # =============================================================================
948 # solver settings:
949 # =============================================================================
950 def setSolverMethod(self,solver=None,preconditioner=None):
951 """
952 sets a new solver
953
954 @param solver: sets a new solver method.
955 @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}
956 @param preconditioner: sets a new solver method.
957 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}, L{GS}
958 """
959 if solver==None: solver=self.__solver_method
960 if preconditioner==None: preconditioner=self.__preconditioner
961 if solver==None: solver=self.DEFAULT
962 if preconditioner==None: preconditioner=self.DEFAULT
963 if not (solver,preconditioner)==self.getSolverMethod():
964 self.__solver_method=solver
965 self.__preconditioner=preconditioner
966 self.__checkMatrixType()
967 self.trace("New solver is %s"%self.getSolverMethodName())
968
969 def getSolverMethodName(self):
970 """
971 returns the name of the solver currently used
972
973 @return: the name of the solver currently used.
974 @rtype: C{string}
975 """
976
977 m=self.getSolverMethod()
978 p=self.getSolverPackage()
979 method=""
980 if m[0]==self.DEFAULT: method="DEFAULT"
981 elif m[0]==self.DIRECT: method= "DIRECT"
982 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
983 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
984 elif m[0]==self.PCG: method= "PCG"
985 elif m[0]==self.TFQMR: method= "TFQMR"
986 elif m[0]==self.MINRES: method= "MINRES"
987 elif m[0]==self.CR: method= "CR"
988 elif m[0]==self.CGS: method= "CGS"
989 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
990 elif m[0]==self.SSOR: method= "SSOR"
991 elif m[0]==self.GMRES: method= "GMRES"
992 elif m[0]==self.PRES20: method= "PRES20"
993 elif m[0]==self.LUMPING: method= "LUMPING"
994 elif m[0]==self.AMG: method= "AMG"
995 if m[1]==self.DEFAULT: method+="+DEFAULT"
996 elif m[1]==self.JACOBI: method+= "+JACOBI"
997 elif m[1]==self.ILU0: method+= "+ILU0"
998 elif m[1]==self.ILUT: method+= "+ILUT"
999 elif m[1]==self.SSOR: method+= "+SSOR"
1000 elif m[1]==self.AMG: method+= "+AMG"
1001 elif m[1]==self.RILU: method+= "+RILU"
1002 elif m[1]==self.GS: method+= "+GS"
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