/[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 425 - (show annotations)
Tue Jan 10 04:10:39 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 110594 byte(s)
The sparse solver can be called by paso now. 

the building has been change to reduce some code redundancy:
now all scons SCscripts are importing scons/esys_options.py which
imports platform specific settings. 



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26