/[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 6799 - (show annotations)
Mon Mar 25 05:53:58 2019 UTC (4 weeks ago) by aellery
File MIME type: text/x-python
File size: 145875 byte(s)
I have rewritten the solverbuddy. Briefly:

1. The remaining AMG code has been removed from PASO.
2. If Trilinos is available, escript will now use it by default.
3. eScript will use a direct solver by default, (if one is available,) when solving 2 dimensional problems and an iterative solver, by default, when solving 3 dimensional problems. This can be changed by a user by manually specifying which solver to use.
4. There is a new option available, setHermitian(), that allows a user to specify when a coefficient matrix is Hermitian.
5. Symmetry information is always passed onto the Trilinos solver when this information is relevant.
6. All tests have been updated, when relevant, to reflect these changes.
7. I fixed a couple of undocumented bugs.


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