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