/[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 2197 - (show annotations)
Thu Jan 8 05:49:16 2009 UTC (10 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 122577 byte(s)
modifications to the way cosntraints are handeled in the Transport problem


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26