/[escript]/branches/diaplayground/escriptcore/py_src/linearPDEs.py
ViewVC logotype

Contents of /branches/diaplayground/escriptcore/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (show annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 140422 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

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