/[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 431 - (show annotations)
Fri Jan 13 05:07:10 2006 UTC (13 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 110838 byte(s)
new implementation of ILU0
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(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1536 if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"
1537 if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"
1538 if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"
1539 mat=self.__getNewOperator()
1540 self.getDomain().addPDEToSystem(mat,escript.Data(), \
1541 self.getCoefficientOfGeneralPDE("A"), \
1542 self.getCoefficientOfGeneralPDE("B"), \
1543 self.getCoefficientOfGeneralPDE("C"), \
1544 self.getCoefficientOfGeneralPDE("D"), \
1545 escript.Data(), \
1546 escript.Data(), \
1547 self.getCoefficientOfGeneralPDE("d"), \
1548 escript.Data(),\
1549 self.getCoefficientOfGeneralPDE("d_contact"), \
1550 escript.Data())
1551 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1552 del mat
1553 self.trace("New lumped operator has been built.")
1554 self.__operator_is_Valid=True
1555 if not self.__righthandside_isValid:
1556 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1557 self.getCoefficientOfGeneralPDE("X"), \
1558 self.getCoefficientOfGeneralPDE("Y"),\
1559 self.getCoefficientOfGeneralPDE("y"),\
1560 self.getCoefficientOfGeneralPDE("y_contact"))
1561 self.trace("New right hand side as been built.")
1562 self.__righthandside_isValid=True
1563 else:
1564 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1565 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1566 self.getCoefficientOfGeneralPDE("A"), \
1567 self.getCoefficientOfGeneralPDE("B"), \
1568 self.getCoefficientOfGeneralPDE("C"), \
1569 self.getCoefficientOfGeneralPDE("D"), \
1570 self.getCoefficientOfGeneralPDE("X"), \
1571 self.getCoefficientOfGeneralPDE("Y"), \
1572 self.getCoefficientOfGeneralPDE("d"), \
1573 self.getCoefficientOfGeneralPDE("y"), \
1574 self.getCoefficientOfGeneralPDE("d_contact"), \
1575 self.getCoefficientOfGeneralPDE("y_contact"))
1576 self.__applyConstraint()
1577 self.__righthandside=self.copyConstraint(self.__righthandside)
1578 self.trace("New system has been built.")
1579 self.__operator_is_Valid=True
1580 self.__righthandside_isValid=True
1581 elif not self.__righthandside_isValid:
1582 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1583 self.getCoefficientOfGeneralPDE("X"), \
1584 self.getCoefficientOfGeneralPDE("Y"),\
1585 self.getCoefficientOfGeneralPDE("y"),\
1586 self.getCoefficientOfGeneralPDE("y_contact"))
1587 self.__righthandside=self.copyConstraint(self.__righthandside)
1588 self.trace("New right hand side has been built.")
1589 self.__righthandside_isValid=True
1590 elif not self.__operator_is_Valid:
1591 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1592 self.getCoefficientOfGeneralPDE("A"), \
1593 self.getCoefficientOfGeneralPDE("B"), \
1594 self.getCoefficientOfGeneralPDE("C"), \
1595 self.getCoefficientOfGeneralPDE("D"), \
1596 escript.Data(), \
1597 escript.Data(), \
1598 self.getCoefficientOfGeneralPDE("d"), \
1599 escript.Data(),\
1600 self.getCoefficientOfGeneralPDE("d_contact"), \
1601 escript.Data())
1602 self.__applyConstraint()
1603 self.trace("New operator has been built.")
1604 self.__operator_is_Valid=True
1605 return (self.__operator,self.__righthandside)
1606
1607
1608 class Poisson(LinearPDE):
1609 """
1610 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1611
1612 M{-grad(grad(u)[j])[j] = f}
1613
1614 with natural boundary conditons
1615
1616 M{n[j]*grad(u)[j] = 0 }
1617
1618 and constraints:
1619
1620 M{u=0} where M{q>0}
1621
1622 """
1623
1624 def __init__(self,domain,debug=False):
1625 """
1626 initializes a new Poisson equation
1627
1628 @param domain: domain of the PDE
1629 @type domain: L{Domain<escript.Domain>}
1630 @param debug: if True debug informations are printed.
1631
1632 """
1633 super(Poisson, self).__init__(domain,1,1,debug)
1634 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1635 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1636 self.setSymmetryOn()
1637
1638 def setValue(self,**coefficients):
1639 """
1640 sets new values to coefficients
1641
1642 @param coefficients: new values assigned to coefficients
1643 @keyword f: value for right hand side M{f}
1644 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1645 @keyword q: mask for location of constraints
1646 @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>}
1647 depending of reduced order is used for the representation of the equation.
1648 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1649 """
1650 super(Poisson, self).setValue(**coefficients)
1651
1652 def getCoefficientOfGeneralPDE(self,name):
1653 """
1654 return the value of the coefficient name of the general PDE
1655 @param name: name of the coefficient requested.
1656 @type name: C{string}
1657 @return: the value of the coefficient name
1658 @rtype: L{Data<escript.Data>}
1659 @raise IllegalCoefficient: if name is not one of coefficients
1660 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}.
1661 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1662 """
1663 if name == "A" :
1664 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1665 elif name == "B" :
1666 return escript.Data()
1667 elif name == "C" :
1668 return escript.Data()
1669 elif name == "D" :
1670 return escript.Data()
1671 elif name == "X" :
1672 return escript.Data()
1673 elif name == "Y" :
1674 return self.getCoefficient("f")
1675 elif name == "d" :
1676 return escript.Data()
1677 elif name == "y" :
1678 return escript.Data()
1679 elif name == "d_contact" :
1680 return escript.Data()
1681 elif name == "y_contact" :
1682 return escript.Data()
1683 elif name == "r" :
1684 return escript.Data()
1685 elif name == "q" :
1686 return self.getCoefficient("q")
1687 else:
1688 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1689
1690 class Helmholtz(LinearPDE):
1691 """
1692 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1693
1694 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1695
1696 with natural boundary conditons
1697
1698 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1699
1700 and constraints:
1701
1702 M{u=r} where M{q>0}
1703
1704 """
1705
1706 def __init__(self,domain,debug=False):
1707 """
1708 initializes a new Poisson equation
1709
1710 @param domain: domain of the PDE
1711 @type domain: L{Domain<escript.Domain>}
1712 @param debug: if True debug informations are printed.
1713
1714 """
1715 super(Helmholtz, self).__init__(domain,1,1,debug)
1716 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1717 "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1718 "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1719 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1720 "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1721 "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1722 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1723 self.setSymmetryOn()
1724
1725 def setValue(self,**coefficients):
1726 """
1727 sets new values to coefficients
1728
1729 @param coefficients: new values assigned to coefficients
1730 @keyword omega: value for coefficient M{S{omega}}
1731 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1732 @keyword k: value for coefficeint M{k}
1733 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1734 @keyword f: value for right hand side M{f}
1735 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1736 @keyword alpha: value for right hand side M{S{alpha}}
1737 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1738 @keyword g: value for right hand side M{g}
1739 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1740 @keyword r: prescribed values M{r} for the solution in constraints.
1741 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1742 depending of reduced order is used for the representation of the equation.
1743 @keyword q: mask for location of constraints
1744 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1745 depending of reduced order is used for the representation of the equation.
1746 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1747 """
1748 super(Helmholtz, self).setValue(**coefficients)
1749
1750 def getCoefficientOfGeneralPDE(self,name):
1751 """
1752 return the value of the coefficient name of the general PDE
1753
1754 @param name: name of the coefficient requested.
1755 @type name: C{string}
1756 @return: the value of the coefficient name
1757 @rtype: L{Data<escript.Data>}
1758 @raise IllegalCoefficient: if name is not one of coefficients
1759 "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}.
1760 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1761 """
1762 if name == "A" :
1763 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
1764 elif name == "B" :
1765 return escript.Data()
1766 elif name == "C" :
1767 return escript.Data()
1768 elif name == "D" :
1769 return self.getCoefficient("omega")
1770 elif name == "X" :
1771 return escript.Data()
1772 elif name == "Y" :
1773 return self.getCoefficient("f")
1774 elif name == "d" :
1775 return self.getCoefficient("alpha")
1776 elif name == "y" :
1777 return self.getCoefficient("g")
1778 elif name == "d_contact" :
1779 return escript.Data()
1780 elif name == "y_contact" :
1781 return escript.Data()
1782 elif name == "r" :
1783 return self.getCoefficient("r")
1784 elif name == "q" :
1785 return self.getCoefficient("q")
1786 else:
1787 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1788
1789 class LameEquation(LinearPDE):
1790 """
1791 Class to define a Lame equation problem:
1792
1793 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] }
1794
1795 with natural boundary conditons:
1796
1797 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] }
1798
1799 and constraints:
1800
1801 M{u[i]=r[i]} where M{q[i]>0}
1802
1803 """
1804
1805 def __init__(self,domain,debug=False):
1806 super(LameEquation, self).__init__(domain,\
1807 domain.getDim(),domain.getDim(),debug)
1808 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1809 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1810 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1811 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1812 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1813 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1814 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1815 self.setSymmetryOn()
1816
1817 def setValue(self,**coefficients):
1818 """
1819 sets new values to coefficients
1820
1821 @param coefficients: new values assigned to coefficients
1822 @keyword lame_mu: value for coefficient M{S{mu}}
1823 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1824 @keyword lame_lambda: value for coefficient M{S{lambda}}
1825 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1826 @keyword F: value for internal force M{F}
1827 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
1828 @keyword sigma: value for initial stress M{S{sigma}}
1829 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
1830 @keyword f: value for extrenal force M{f}
1831 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
1832 @keyword r: prescribed values M{r} for the solution in constraints.
1833 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1834 depending of reduced order is used for the representation of the equation.
1835 @keyword q: mask for location of constraints
1836 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1837 depending of reduced order is used for the representation of the equation.
1838 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1839 """
1840 super(LameEquation, self).setValue(**coefficients)
1841
1842 def getCoefficientOfGeneralPDE(self,name):
1843 """
1844 return the value of the coefficient name of the general PDE
1845
1846 @param name: name of the coefficient requested.
1847 @type name: C{string}
1848 @return: the value of the coefficient name
1849 @rtype: L{Data<escript.Data>}
1850 @raise IllegalCoefficient: if name is not one of coefficients
1851 "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}.
1852 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1853 """
1854 if name == "A" :
1855 out =self.createCoefficientOfGeneralPDE("A")
1856 for i in range(self.getDim()):
1857 for j in range(self.getDim()):
1858 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1859 out[i,j,j,i] += self.getCoefficient("lame_mu")
1860 out[i,j,i,j] += self.getCoefficient("lame_mu")
1861 return out
1862 elif name == "B" :
1863 return escript.Data()
1864 elif name == "C" :
1865 return escript.Data()
1866 elif name == "D" :
1867 return escript.Data()
1868 elif name == "X" :
1869 return self.getCoefficient("sigma")
1870 elif name == "Y" :
1871 return self.getCoefficient("F")
1872 elif name == "d" :
1873 return escript.Data()
1874 elif name == "y" :
1875 return self.getCoefficient("f")
1876 elif name == "d_contact" :
1877 return escript.Data()
1878 elif name == "y_contact" :
1879 return escript.Data()
1880 elif name == "r" :
1881 return self.getCoefficient("r")
1882 elif name == "q" :
1883 return self.getCoefficient("q")
1884 else:
1885 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1886
1887 class AdvectivePDE(LinearPDE):
1888 """
1889 In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}
1890 up-winding has been used. The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.
1891
1892 In the following we set
1893
1894 M{Z[j]=C[j]-B[j]}
1895
1896 or
1897
1898 M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1899
1900 To measure the dominance of the advective terms over the diffusive term M{A} the
1901 X{Pelclet number} M{P} is used. It is defined as
1902
1903 M{P=h|Z|/(2|A|)}
1904
1905 where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1906 from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1907
1908 From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
1909
1910 M{S{Xi}=S{xi}(P) h/|Z|}
1911
1912 where M{S{xi}} is a suitable function of the Peclet number.
1913
1914 In the case of a single PDE the coefficient are up-dated in the following way:
1915 - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}
1916 - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}
1917 - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}
1918 - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}
1919
1920 Similar for the case of a systems of PDEs:
1921 - 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]}
1922 - M{B[i,j,k] S{<-} B[i,j,k] + S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}
1923 - M{C[i,k,l] S{<-} C[i,k,l] + S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}
1924 - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi} * Y[p] * Z[m,i,j]}
1925
1926 where M{S{delta}} is L{kronecker}.
1927 Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}
1928 but with the intension to stabilize the solution.
1929
1930 """
1931 def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1932 """
1933 creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1934
1935 @param domain: domain of the PDE
1936 @type domain: L{Domain<escript.Domain>}
1937 @param numEquations: number of equations. If numEquations==None the number of equations
1938 is exracted from the PDE coefficients.
1939 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1940 is exracted from the PDE coefficients.
1941 @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1942 M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1943 @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1944 @param debug: if True debug informations are printed.
1945 """
1946 super(AdvectivePDE, self).__init__(domain,\
1947 numEquations,numSolutions,debug)
1948 if xi==None:
1949 self.__xi=AdvectivePDE.ELMAN_RAMAGE
1950 else:
1951 self.__xi=xi
1952 self.__Xi=escript.Data()
1953
1954 def setValue(**coefficients):
1955 """
1956 sets new values to coefficients
1957
1958 @param coefficients: new values assigned to coefficients
1959 @keyword A: value for coefficient A.
1960 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1961 @keyword B: value for coefficient B
1962 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1963 @keyword C: value for coefficient C
1964 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1965 @keyword D: value for coefficient D
1966 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1967 @keyword X: value for coefficient X
1968 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1969 @keyword Y: value for coefficient Y
1970 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1971 @keyword d: value for coefficient d
1972 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1973 @keyword y: value for coefficient y
1974 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1975 @keyword d_contact: value for coefficient d_contact
1976 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1977 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1978 @keyword y_contact: value for coefficient y_contact
1979 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1980 or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1981 @keyword r: values prescribed to the solution at the locations of constraints
1982 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1983 depending of reduced order is used for the solution.
1984 @keyword q: mask for location of constraints
1985 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1986 depending of reduced order is used for the representation of the equation.
1987 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1988
1989 """
1990 if "A" in coefficients.keys() or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
1991 super(AdvectivePDE, self).setValue(**coefficients)
1992
1993 def ELMAN_RAMAGE(self,P):
1994 """
1995 Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1996 This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
1997 - M{S{xi}(P)=0} for M{P<1}
1998 - M{S{xi}(P)=(1-1/P)/2} otherwise
1999
2000 @param P: Preclet number
2001 @type P: L{Scalar<escript.Scalar>}
2002 @return: up-wind weightimg factor
2003 @rtype: L{Scalar<escript.Scalar>}
2004 """
2005 return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
2006
2007 def SIMPLIFIED_BROOK_HUGHES(self,P):
2008 """
2009 Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2010 The original methods is
2011
2012 M{S{xi}(P)=coth(P)-1/P}
2013
2014 As the evaluation of M{coth} is expensive we are using the approximation:
2015
2016 - M{S{xi}(P)=P/3} where M{P<3}
2017 - M{S{xi}(P)=1/2} otherwise
2018
2019 @param P: Preclet number
2020 @type P: L{Scalar<escript.Scalar>}
2021 @return: up-wind weightimg factor
2022 @rtype: L{Scalar<escript.Scalar>}
2023 """
2024 c=util.whereNegative(P-3.)
2025 return P/6.*c+1./2.*(1.-c)
2026
2027 def HALF(self,P):
2028 """
2029 Predefined function to set value M{1/2} for M{S{xi}}
2030
2031 @param P: Preclet number
2032 @type P: L{Scalar<escript.Scalar>}
2033 @return: up-wind weightimg factor
2034 @rtype: L{Scalar<escript.Scalar>}
2035 """
2036 return escript.Scalar(0.5,P.getFunctionSpace())
2037
2038 def __calculateXi(self,peclet_factor,flux,h):
2039 flux=util.Lsup(flux)
2040 if flux_max>0.:
2041 return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2042 else:
2043 return 0.
2044
2045 def __getXi(self):
2046 if self.__Xi.isEmpty():
2047 B=self.getCoefficient("B")
2048 C=self.getCoefficient("C")
2049 A=self.getCoefficient("A")
2050 h=self.getDomain().getSize()
2051 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2052 if not C.isEmpty() or not B.isEmpty():
2053 if not C.isEmpty() and not B.isEmpty():
2054 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2055 if self.getNumEquations()>1:
2056 if self.getNumSolutions()>1:
2057 for i in range(self.getNumEquations()):
2058 for k in range(self.getNumSolutions()):
2059 for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2060 # flux=C-util.reorderComponents(B,[0,2,1])
2061 else:
2062 for i in range(self.getNumEquations()):
2063 for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2064 # flux=C-B
2065 else:
2066 if self.getNumSolutions()>1:
2067 for k in range(self.getNumSolutions()):
2068 for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2069 # flux=C-util.reorderComponents(B,[1,0])
2070 else:
2071 for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2072 #flux=C-B
2073 length_of_flux=util.sqrt(flux2)
2074 elif C.isEmpty():
2075 length_of_flux=util.length(B)
2076 #flux=B
2077 else:
2078 length_of_flux=util.length(C)
2079 #flux=C
2080
2081 #length_of_flux=util.length(flux)
2082 flux_max=util.Lsup(length_of_flux)
2083 if flux_max>0.:
2084 # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2085 length_of_A=util.length(A)
2086 A_max=util.Lsup(length_of_A)
2087 if A_max>0:
2088 inv_A=1./(length_of_A+A_max*self.__TOL)
2089 else:
2090 inv_A=1./self.__TOL
2091 peclet_number=length_of_flux*h/2*inv_A
2092 xi=self.__xi(peclet_number)
2093 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2094 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2095 return self.__Xi
2096
2097
2098 def getCoefficientOfGeneralPDE(self,name):
2099 """
2100 return the value of the coefficient name of the general PDE
2101
2102 @param name: name of the coefficient requested.
2103 @type name: C{string}
2104 @return: the value of the coefficient name
2105 @rtype: L{Data<escript.Data>}
2106 @raise IllegalCoefficient: if name is not one of coefficients
2107 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}.
2108 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2109 """
2110 if not self.getNumEquations() == self.getNumSolutions():
2111 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
2112
2113 if name == "A" :
2114 A=self.getCoefficient("A")
2115 B=self.getCoefficient("B")
2116 C=self.getCoefficient("C")
2117 if B.isEmpty() and C.isEmpty():
2118 Aout=A
2119 else:
2120 if A.isEmpty():
2121 Aout=self.createNewCoefficient("A")
2122 else:
2123 Aout=A[:]
2124 Xi=self.__getXi()
2125 if self.getNumEquations()>1:
2126 for i in range(self.getNumEquations()):
2127 for j in range(self.getDim()):
2128 for k in range(self.getNumSolutions()):
2129 for l in range(self.getDim()):
2130 if not C.isEmpty() and not B.isEmpty():
2131 # tmp=C-util.reorderComponents(B,[0,2,1])
2132 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2133 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])
2134 elif C.isEmpty():
2135 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2136 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2137 else:
2138 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2139 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2140 else:
2141 for j in range(self.getDim()):
2142 for l in range(self.getDim()):
2143 if not C.isEmpty() and not B.isEmpty():
2144 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
2145 elif C.isEmpty():
2146 Aout[j,l]+=Xi*B[j]*B[l]
2147 else:
2148 Aout[j,l]+=Xi*C[j]*C[l]
2149 # if not C.isEmpty() and not B.isEmpty():
2150 # tmp=C-B
2151 # Aout=Aout+Xi*util.outer(tmp,tmp)
2152 # elif C.isEmpty():
2153 # Aout=Aout+Xi*util.outer(B,B)
2154 # else:
2155 # Aout=Aout+Xi*util.outer(C,C)
2156 return Aout
2157 elif name == "B" :
2158 B=self.getCoefficient("B")
2159 C=self.getCoefficient("C")
2160 D=self.getCoefficient("D")
2161 if C.isEmpty() or D.isEmpty():
2162 Bout=B
2163 else:
2164 Xi=self.__getXi()
2165 if B.isEmpty():
2166 Bout=self.createNewCoefficient("B")
2167 else:
2168 Bout=B[:]
2169 if self.getNumEquations()>1:
2170 for k in range(self.getNumSolutions()):
2171 for p in range(self.getNumEquations()):
2172 tmp=Xi*D[p,k]
2173 for i in range(self.getNumEquations()):
2174 for j in range(self.getDim()):
2175 Bout[i,j,k]+=tmp*C[p,i,j]
2176 # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2177 else:
2178 tmp=Xi*D
2179 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2180 # Bout=Bout+Xi*D*C
2181 return Bout
2182 elif name == "C" :
2183 B=self.getCoefficient("B")
2184 C=self.getCoefficient("C")
2185 D=self.getCoefficient("D")
2186 if B.isEmpty() or D.isEmpty():
2187 Cout=C
2188 else:
2189 Xi=self.__getXi()
2190 if C.isEmpty():
2191 Cout=self.createNewCoefficient("C")
2192 else:
2193 Cout=C[:]
2194 if self.getNumEquations()>1:
2195 for k in range(self.getNumSolutions()):
2196 for p in range(self.getNumEquations()):
2197 tmp=Xi*D[p,k]
2198 for i in range(self.getNumEquations()):
2199 for l in range(self.getDim()):
2200 Cout[i,k,l]+=tmp*B[p,l,i]
2201 # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2202 else:
2203 tmp=Xi*D
2204 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2205 # Cout=Cout+tmp*D*B
2206 return Cout
2207 elif name == "D" :
2208 return self.getCoefficient("D")
2209 elif name == "X" :
2210 X=self.getCoefficient("X")
2211 Y=self.getCoefficient("Y")
2212 B=self.getCoefficient("B")
2213 C=self.getCoefficient("C")
2214 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
2215 Xout=X
2216 else:
2217 if X.isEmpty():
2218 Xout=self.createNewCoefficient("X")
2219 else:
2220 Xout=X[:]
2221 Xi=self.__getXi()
2222 if self.getNumEquations()>1:
2223 for p in range(self.getNumEquations()):
2224 tmp=Xi*Y[p]
2225 for i in range(self.getNumEquations()):
2226 for j in range(self.getDim()):
2227 if not C.isEmpty() and not B.isEmpty():
2228 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2229 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2230 elif C.isEmpty():
2231 Xout[i,j]-=tmp*B[p,j,i]
2232 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2233 else:
2234 Xout[i,j]+=tmp*C[p,i,j]
2235 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2236 else:
2237 tmp=Xi*Y
2238 for j in range(self.getDim()):
2239 if not C.isEmpty() and not B.isEmpty():
2240 Xout[j]+=tmp*(C[j]-B[j])
2241 # Xout=Xout+Xi*Y*(C-B)
2242 elif C.isEmpty():
2243 Xout[j]-=tmp*B[j]
2244 # Xout=Xout-Xi*Y*B
2245 else:
2246 Xout[j]+=tmp*C[j]
2247 # Xout=Xout+Xi*Y*C
2248 return Xout
2249 elif name == "Y" :
2250 return self.getCoefficient("Y")
2251 elif name == "d" :
2252 return self.getCoefficient("d")
2253 elif name == "y" :
2254 return self.getCoefficient("y")
2255 elif name == "d_contact" :
2256 return self.getCoefficient("d_contact")
2257 elif name == "y_contact" :
2258 return self.getCoefficient("y_contact")
2259 elif name == "r" :
2260 return self.getCoefficient("r")
2261 elif name == "q" :
2262 return self.getCoefficient("q")
2263 else:
2264 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2265
2266 class AdvectionDiffusion(LinearPDE):
2267 """
2268 Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2269
2270 M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}
2271
2272 with natural boundary conditons
2273
2274 M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2275
2276 and constraints:
2277
2278 M{u=r} where M{q>0}
2279
2280 and
2281
2282 M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}
2283
2284 """
2285
2286 def __init__(self,domain,debug=False):
2287 """
2288 initializes a new Poisson equation
2289
2290 @param domain: domain of the PDE
2291 @type domain: L{Domain<escript.Domain>}
2292 @param debug: if True debug informations are printed.
2293
2294 """
2295 super(AdvectionDiffusion, self).__init__(domain,1,1,debug)
2296 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2297 "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2298 "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2299 "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2300 "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2301 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2302 "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2303 "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2304 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2305
2306 def setValue(self,**coefficients):
2307 """
2308 sets new values to coefficients
2309
2310 @param coefficients: new values assigned to coefficients
2311 @keyword omega: value for coefficient M{S{omega}}
2312 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2313 @keyword k: value for coefficient M{k}
2314 @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2315 @keyword v: value for coefficient M{v}
2316 @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2317 @keyword upwind: value for upwind term M{upwind}
2318 @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2319 @keyword f: value for right hand side M{f}
2320 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2321 @keyword alpha: value for right hand side M{S{alpha}}
2322 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2323 @keyword g: value for right hand side M{g}
2324 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2325 @keyword r: prescribed values M{r} for the solution in constraints.
2326 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2327 depending of reduced order is used for the representation of the equation.
2328 @keyword q: mask for location of constraints
2329 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2330 depending of reduced order is used for the representation of the equation.
2331 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2332 """
2333 super(AdvectionDiffusion, self).setValue(**coefficients)
2334
2335 def getCoefficientOfGeneralPDE(self,name):
2336 """
2337 return the value of the coefficient name of the general PDE
2338
2339 @param name: name of the coefficient requested.
2340 @type name: C{string}
2341 @return: the value of the coefficient name
2342 @rtype: L{Data<escript.Data>}
2343 @raise IllegalCoefficient: if name is not one of coefficients
2344 "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}.
2345 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2346 """
2347 if name == "A" :
2348 return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))
2349 elif name == "B" :
2350 return escript.Data()
2351 elif name == "C" :
2352 return self.getCoefficient("v")
2353 elif name == "D" :
2354 return self.getCoefficient("omega")
2355 elif name == "X" :
2356 return escript.Data()
2357 elif name == "Y" :
2358 return self.getCoefficient("f")
2359 elif name == "d" :
2360 return self.getCoefficient("alpha")
2361 elif name == "y" :
2362 return self.getCoefficient("g")
2363 elif name == "d_contact" :
2364 return escript.Data()
2365 elif name == "y_contact" :
2366 return escript.Data()
2367 elif name == "r" :
2368 return self.getCoefficient("r")
2369 elif name == "q" :
2370 return self.getCoefficient("q")
2371 else:
2372 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2373
2374
2375 # $Log$
2376 # Revision 1.14 2005/09/22 01:54:57 jgs
2377 # Merge of development branch dev-02 back to main trunk on 2005-09-22
2378 #
2379 # Revision 1.13 2005/09/15 03:44:19 jgs
2380 # Merge of development branch dev-02 back to main trunk on 2005-09-15
2381 #
2382 # Revision 1.12 2005/09/01 03:31:28 jgs
2383 # Merge of development branch dev-02 back to main trunk on 2005-09-01
2384 #
2385 # Revision 1.11 2005/08/23 01:24:28 jgs
2386 # Merge of development branch dev-02 back to main trunk on 2005-08-23
2387 #
2388 # Revision 1.10 2005/08/12 01:45:36 jgs
2389 # erge of development branch dev-02 back to main trunk on 2005-08-12
2390 #
2391 # Revision 1.9.2.17 2005/09/21 07:03:33 matt
2392 # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2393 # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2394 # modified to instead use portable/cooperative "super" calls to extend base
2395 # class methods.
2396 #
2397 # Revision 1.9.2.16 2005/09/16 01:54:37 matt
2398 # Removed redundant if-loop.
2399 #
2400 # Revision 1.9.2.15 2005/09/14 08:09:18 matt
2401 # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2402 #
2403 # Revision 1.9.2.14 2005/09/07 06:26:16 gross
2404 # the solver from finley are put into the standalone package paso now
2405 #
2406 # Revision 1.9.2.13 2005/08/31 08:45:03 gross
2407 # in the case of lumping no new system is allocated if the constraint is changed.
2408 #
2409 # Revision 1.9.2.12 2005/08/31 07:10:23 gross
2410 # test for Lumping added
2411 #
2412 # Revision 1.9.2.11 2005/08/30 01:53:45 gross
2413 # bug in format fixed.
2414 #
2415 # Revision 1.9.2.10 2005/08/26 07:14:17 gross
2416 # a few more bugs in linearPDE fixed. remaining problem are finley problems
2417 #
2418 # Revision 1.9.2.9 2005/08/26 06:30:45 gross
2419 # fix for reported bug 0000004. test_linearPDE passes a few more tests
2420 #
2421 # Revision 1.9.2.8 2005/08/26 04:30:13 gross
2422 # gneric unit testing for linearPDE
2423 #
2424 # Revision 1.9.2.7 2005/08/25 07:06:50 gross
2425 # linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so
2426 #
2427 # Revision 1.9.2.6 2005/08/24 05:01:24 gross
2428 # problem with resetting the matrix in case of resetting its values to 0 fixed.
2429 #
2430 # Revision 1.9.2.5 2005/08/24 02:03:28 gross
2431 # epydoc mark up partially fixed
2432 #
2433 # Revision 1.9.2.4 2005/08/22 07:11:09 gross
2434 # some problems with LinearPDEs fixed.
2435 #
2436 # Revision 1.9.2.3 2005/08/18 04:48:48 gross
2437 # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
2438 #
2439 # Revision 1.9.2.2 2005/08/18 04:39:32 gross
2440 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
2441 #
2442 # Revision 1.9.2.1 2005/07/29 07:10:27 gross
2443 # new functions in util and a new pde type in linearPDEs
2444 #
2445 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
2446 # Lame equation: (linear elastic, isotropic) added
2447 #
2448 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
2449 # some extensions to modellib and linearPDEs
2450 #
2451 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
2452 # Fixed up some docstrings. Moved module-level functions to top of file so
2453 # that epydoc and doxygen can pick them up properly.
2454 #
2455 # Revision 1.1.2.22 2005/05/12 11:41:30 gross
2456 # some basic Models have been added
2457 #
2458 # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
2459 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
2460 # file so that the AdvectivePDE class is picked up by doxygen. Some
2461 # reformatting of docstrings. Addition of code to make equations come out
2462 # as proper LaTeX.
2463 #
2464 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
2465 # some problems with functionspace and linearPDEs fixed.
2466 #
2467 # Revision 1.1.2.19 2005/03/04 05:27:07 gross
2468 # bug in SystemPattern fixed.
2469 #
2470 # Revision 1.1.2.18 2005/02/08 06:16:45 gross
2471 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
2472 #
2473 # Revision 1.1.2.17 2005/02/08 05:56:19 gross
2474 # Reference Number handling added
2475 #
2476 # Revision 1.1.2.16 2005/02/07 04:41:28 gross
2477 # some function exposed to python to make mesh merging running
2478 #
2479 # Revision 1.1.2.15 2005/02/03 00:14:44 gross
2480 # timeseries add and ESySParameter.py renames esysXML.py for consistence
2481 #
2482 # Revision 1.1.2.14 2005/02/01 06:44:10 gross
2483 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
2484 #
2485 # Revision 1.1.2.13 2005/01/25 00:47:07 gross
2486 # updates in the documentation
2487 #
2488 # Revision 1.1.2.12 2005/01/12 01:28:04 matt
2489 # Added createCoefficient method for linearPDEs.
2490 #
2491 # Revision 1.1.2.11 2005/01/11 01:55:34 gross
2492 # a problem in linearPDE class fixed
2493 #
2494 # Revision 1.1.2.10 2005/01/07 01:13:29 gross
2495 # some bugs in linearPDE fixed
2496 #
2497 # Revision 1.1.2.9 2005/01/06 06:24:58 gross
2498 # some bugs in slicing fixed
2499 #
2500 # Revision 1.1.2.8 2005/01/05 04:21:40 gross
2501 # FunctionSpace checking/matchig in slicing added
2502 #
2503 # Revision 1.1.2.7 2004/12/29 10:03:41 gross
2504 # bug in setValue fixed
2505 #
2506 # Revision 1.1.2.6 2004/12/29 05:29:59 gross
2507 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
2508 #
2509 # Revision 1.1.2.5 2004/12/29 00:18:41 gross
2510 # AdvectivePDE added
2511 #
2512 # Revision 1.1.2.4 2004/12/24 06:05:41 gross
2513 # some changes in linearPDEs to add AdevectivePDE
2514 #
2515 # Revision 1.1.2.3 2004/12/16 00:12:34 gross
2516 # __init__ of LinearPDE does not accept any coefficient anymore
2517 #
2518 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
2519 # *** empty log message ***
2520 #
2521 # Revision 1.1.2.1 2004/12/12 22:53:47 gross
2522 # linearPDE has been renamed LinearPDE
2523 #
2524 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
2525 # GMRES added
2526 #
2527 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
2528 # options for GMRES and PRES20 added
2529 #
2530 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
2531 # some small changes
2532 #
2533 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
2534 # Finley solves 4M unknowns now
2535 #
2536 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
2537 # poisson solver added
2538 #
2539 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
2540 # 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
2541 #
2542 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
2543 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
2544 #
2545 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
2546 # initial import of project esys2
2547 #
2548 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
2549 # committing Lutz's and Paul's changes to brach jgs
2550 #
2551 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
2552 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
2553 #
2554 # Revision 1.3 2004/09/23 00:53:23 jgs
2555 # minor fixes
2556 #
2557 # Revision 1.1 2004/08/28 12:58:06 gross
2558 # SimpleSolve is not running yet: problem with == of functionsspace
2559 #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26