/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3051 - (show annotations)
Tue Jun 29 00:45:49 2010 UTC (9 years, 4 months ago) by lgao
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 152040 byte(s)
Hybrid MPI/OpenMP versioned Gauss_Seidel preconditioner is added. 
To use it, use "SolverOptions.GAUSS_SEIDEL_MPI" in python script. 


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