/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years, 1 month ago) by ksteube
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 112970 byte(s)
Copyright updated in all python files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26