/[escript]/trunk/escriptcore/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk/escriptcore/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6670 - (show annotations)
Mon Apr 30 03:03:17 2018 UTC (3 years, 1 month ago) by uqagarro
File MIME type: text/x-python
File size: 144523 byte(s)
Raising exception for worng q type on linearPDE call

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