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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26