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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26