/[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 3517 - (show annotations)
Fri May 20 01:16:41 2011 UTC (8 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 163852 byte(s)
Merged up to and including revision 3514 from trunk and implemented
symbolic inverse.

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