/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


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