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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (hide annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 167699 byte(s)
Fast forward to latest trunk revision 3788.

1 gross 3094 # -*- coding: utf-8 -*-
2 ksteube 1809
3     ########################################################
4 ksteube 1312 #
5 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
6 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
7     # http://www.uq.edu.au/esscc
8 ksteube 1312 #
9 ksteube 1809 # 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 ksteube 1312 #
13 ksteube 1809 ########################################################
14 ksteube 1312
15 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
22 ksteube 1809
23 jgs 149 """
24 jgs 151 The module provides an interface to define and solve linear partial
25 jfenwick 2625 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 caltinay 2169 by its advective terms.
31 jgs 102
32 jfenwick 2625 :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 jgs 102 """
39    
40 gross 1417 import math
41 caltinay 3789 from . import escript
42     from . import util
43 jfenwick 2455 import numpy
44 jgs 102
45 jgs 149 __author__="Lutz Gross, l.gross@uq.edu.au"
46 jgs 102
47 jgs 149
48 gross 2470 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 jfenwick 2625 ::
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 gross 2470
64 jfenwick 2625 :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 gross 3071 :cvar TFQMR: Transpose Free Quasi Minimal Residual method
72 jfenwick 2625 :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 gross 3379 :cvar ROWSUM_LUMPING: Matrix lumping using row sum
79     :cvar HRZ_LUMPING: Matrix lumping using the HRZ approach
80 jfenwick 2625 :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 artak 2760 :cvar MKL: Intel's MKL solver library
86 jfenwick 2625 :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 artak 2760 :cvar AMLI: Algebraic Multi Level Iteration
91 jfenwick 2625 :cvar REC_ILU: recursive ILU0
92     :cvar RILU: relaxed ILU0
93 gross 3094 :cvar GAUSS_SEIDEL: Gauss-Seidel preconditioner
94 jfenwick 2625 :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 artak 2816 :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 caltinay 3517 :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
101 jfenwick 2625 :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 gross 3440 :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 caltinay 3517 :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 gross 2470 """
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 gross 3379 ROWSUM_LUMPING= 13
128     HRZ_LUMPING= 14
129 gross 2470 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 artak 2524 MIN_COARSE_MATRIX_SIZE=37
152 artak 2816 AMLI=38
153     STANDARD_COARSENING=39
154 gross 3440 CLASSIC_INTERPOLATION_WITH_FF_COUPLING=50
155     CLASSIC_INTERPOLATION=51
156     DIRECT_INTERPOLATION=52
157 caltinay 3517 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 artak 2760
164 gross 2470 def __init__(self):
165     self.setLevelMax()
166     self.setCoarseningThreshold()
167 artak 2828 self.setSmoother()
168 gross 2470 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 artak 2524 self.setMinCoarseMatrixSize()
190 gross 2470 self.setRelaxationFactor()
191 gross 3094 self.setLocalPreconditionerOff()
192 gross 2470 self.resetDiagnostics(all=True)
193 gross 3283 self.setMinCoarseMatrixSparsity()
194     self.setNumRefinements()
195     self.setNumCoarseMatrixRefinements()
196 gross 3352 self.setUsePanel()
197 gross 3402 self.setDiagonalDominanceThreshold()
198 gross 3440 self.setAMGInterpolation()
199 caltinay 3517 self.setCycleType()
200 gross 3283
201 gross 2470
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 gross 3094 # out+="\nInner tolerance = %e"%self.getInnerTolerance()
216     # out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
217 gross 2474
218 gross 2470 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 gross 3283 out+="\nMaximum number of levels = %s"%self.getLevelMax()
226 gross 2470 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
227 gross 3094 out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
228 gross 2470 out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
229 gross 3094 out+="\nApply preconditioner locally = %s"%self.useLocalPreconditioner()
230 gross 2470 if self.getPreconditioner() == self.AMG:
231 gross 3283 out+="\nMaximum number of levels = %s"%self.getLevelMax()
232 gross 3352 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
233     out+="\nMinimal sparsity on coarsest level = %e"%(self.getMinCoarseMatrixSparsity(),)
234 gross 3283 out+="\nSmoother = %s"%self.getName(self.getSmoother())
235 artak 2524 out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
236 gross 2470 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
237 gross 3283 out+="\nNumber of refinement steps in coarsest level solver = %d"%(self.getNumCoarseMatrixRefinements(),)
238 gross 3352 out+="\nUse node panel = %s"%(self.usePanel())
239 gross 3440 out+="\nInterpolation = %s"%(self.getName(self.getAMGInterpolation()))
240 gross 3402 out+="\nThreshold for diagonal dominant rows = %s"%(setDiagonalDominanceThreshold(),)
241 gross 3352
242 artak 2760 if self.getPreconditioner() == self.AMLI:
243 gross 3283 out+="\nMaximum number of levels = %s"%self.getLevelMax()
244 artak 2760 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 caltinay 3517 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 artak 2760 if self.getPreconditioner() == self.GAUSS_SEIDEL:
257 gross 2470 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 jfenwick 2625 :param key: a valid key
270 gross 2470 """
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 gross 3379 if key == self.ROWSUM_LUMPING: return "ROWSUM_LUMPING"
284     if key == self.HRZ_LUMPING: return "HRZ_LUMPING"
285 gross 2470 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 artak 2760 if key == self.AMLI: return "AMLI"
294 gross 2470 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 artak 2816 if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
307 gross 2470 if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
308     if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
309 gross 3440 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 caltinay 3517 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 gross 2470
319     def resetDiagnostics(self,all=False):
320     """
321     resets the diagnostics
322    
323 jfenwick 2625 :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
324     :type all: ``bool``
325 gross 2470 """
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 lgao 2596 self.__net_time=None
332 gross 2470 self.__residual_norm=None
333     self.__converged=None
334 gross 3094 self.__preconditioner_size=-1
335 gross 3111 self.__time_step_backtracking_used=None
336 gross 2470 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 lgao 2596 self.__cum_net_time=0
342 gross 2470
343     def _updateDiagnostics(self, name, value):
344     """
345     Updates diagnostic information
346    
347 jfenwick 2625 :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 gross 2470 """
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 lgao 2596 self.__num_level=int(value)
357 gross 2470 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 lgao 2596 if name == "net_time":
367     self.__net_time=float(value)
368     self.__cum_net_time+=self.__net_time
369 gross 2470 if name == "residual_norm":
370     self.__residual_norm=float(value)
371     if name == "converged":
372     self.__converged = (value == True)
373 gross 3111 if name == "time_step_backtracking_used":
374     self.__time_step_backtracking_used = (value == True)
375 gross 3402 if name == "coarse_level_sparsity":
376     self.__coarse_level_sparsity = value
377     if name == "num_coarse_unknowns":
378     self.__num_coarse_unknowns= value
379 gross 2470 def getDiagnostics(self, name):
380     """
381 jfenwick 2625 Returns the diagnostic information ``name``. Possible values are:
382    
383     - "num_iter": the number of iteration steps
384 gross 2470 - "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 lgao 2596 - "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 gross 3094 - "preconditioner_size": size of preconditioner [Bytes]
395 gross 3402 - "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 jfenwick 2625
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 gross 2470 """
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 lgao 2596 elif name == "net_time": return self.__net_time
416     elif name == "cum_net_time": return self.__cum_net_time
417 gross 2470 elif name == "residual_norm": return self.__residual_norm
418     elif name == "converged": return self.__converged
419 gross 3094 elif name == "preconditioner_size": return self.__preconditioner_size
420 gross 3111 elif name == "time_step_backtracking_used": return self.__time_step_backtracking_used
421 gross 3402 elif name == "coarse_level_sparsity": return self.__coarse_level_sparsity
422     elif name == "num_coarse_unknowns": return self.__num_coarse_unknowns
423 gross 2470 else:
424 caltinay 3789 raise ValueError("unknown diagnostic item %s"%name)
425 gross 2470 def hasConverged(self):
426     """
427 jfenwick 2625 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 gross 2470 """
430     return self.getDiagnostics("converged")
431     def setCoarsening(self,method=0):
432     """
433 caltinay 3517 Sets the key of the coarsening method to be applied in AMG or AMLI or BoomerAMG
434 gross 2470
435 jfenwick 2625 :param method: selects the coarsening method .
436 caltinay 3517 :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 gross 2470 """
438 gross 2474 if method==None: method=0
439 caltinay 3517 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 caltinay 3789 raise ValueError("unknown coarsening method %s"%method)
441 gross 2470 self.__coarsening=method
442 artak 2524
443 gross 2470 def getCoarsening(self):
444     """
445 caltinay 3517 Returns the key of the coarsening algorithm to be applied AMG or AMLI or BoomerAMG
446 gross 2470
447 caltinay 3517 :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 gross 2470 """
449     return self.__coarsening
450 artak 2524
451 gross 3433 def setMinCoarseMatrixSize(self,size=None):
452 artak 2524 """
453 gross 3094 Sets the minumum size of the coarsest level matrix in AMG or AMLI
454 artak 2524
455 jfenwick 2625 :param size: minumum size of the coarsest level matrix .
456     :type size: positive ``int`` or ``None``
457 artak 2524 """
458 gross 3094 if size==None: size=500
459 artak 2524 size=int(size)
460     if size<0:
461 caltinay 3789 raise ValueError("minumum size of the coarsest level matrix must be non-negative.")
462 artak 2524 self.__MinCoarseMatrixSize=size
463    
464     def getMinCoarseMatrixSize(self):
465     """
466 gross 3094 Returns the minumum size of the coarsest level matrix in AMG or AMLI
467 artak 2524
468 jfenwick 2625 :rtype: ``int``
469 artak 2524 """
470     return self.__MinCoarseMatrixSize
471    
472 gross 2470 def setPreconditioner(self, preconditioner=10):
473     """
474     Sets the preconditioner to be used.
475    
476 jfenwick 2625 :param preconditioner: key of the preconditioner to be used.
477 gross 3094 :type preconditioner: in `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
478 caltinay 3517 `SolverOptions.AMG`, `SolverOptions.AMLI`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`, `SolverOptions.BOOMERAMG`,
479 jfenwick 2625 `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 gross 2470 """
482 gross 2474 if preconditioner==None: preconditioner=10
483 gross 3094 if not preconditioner in [ SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
484 caltinay 3517 SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU, SolverOptions.BOOMERAMG,
485 gross 2470 SolverOptions.NO_PRECONDITIONER] :
486 caltinay 3789 raise ValueError("unknown preconditioner %s"%preconditioner)
487 gross 2470 self.__preconditioner=preconditioner
488     def getPreconditioner(self):
489     """
490     Returns key of the preconditioner to be used.
491    
492 gross 3094 :rtype: in the list `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`, SolverOptions.AMLI,
493 caltinay 3517 `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`, `Solveroptions.BOOMERAMG`,
494 jfenwick 2625 `SolverOptions.NO_PRECONDITIONER`
495 gross 2470 """
496     return self.__preconditioner
497 gross 3094 def setSmoother(self, smoother=None):
498 artak 2828 """
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 artak 2835 if smoother==None: smoother=28
506 artak 2828 if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
507 caltinay 3789 raise ValueError("unknown smoother %s"%smoother)
508 artak 2828 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 gross 3440
517 gross 2470 def setSolverMethod(self, method=0):
518     """
519 jfenwick 2625 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 gross 2470 solver should be used.
522    
523 jfenwick 2625 :param method: key of the solver method to be used.
524     :type method: in `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
525 gross 3094 `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
526 gross 3379 `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
527 gross 3094 `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
528 jfenwick 2625 :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 gross 2470 """
530 gross 2474 if method==None: method=0
531 gross 2470 if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
532 gross 3094 SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,
533 gross 3379 SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
534     SolverOptions.ITERATIVE,
535 gross 3094 SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:
536 caltinay 3789 raise ValueError("unknown solver method %s"%method)
537 gross 2470 self.__method=method
538     def getSolverMethod(self):
539     """
540     Returns key of the solver method to be used.
541    
542 jfenwick 2625 :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
543 gross 3094 `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
544 gross 3379 `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
545 gross 3094 `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
546 gross 2470 """
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 jfenwick 2625 :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 gross 2470 """
557 gross 2474 if package==None: package=0
558 gross 2470 if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
559 caltinay 3789 raise ValueError("unknown solver package %s"%package)
560 gross 2470 self.__package=package
561     def getPackage(self):
562     """
563     Returns the solver package key
564    
565 jfenwick 2625 :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
566 gross 2470 """
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 jfenwick 2625 :param ordering: selects the reordering strategy.
574 gross 3094 :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
575 gross 2470 """
576     if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
577 caltinay 3789 raise ValueError("unknown reordering strategy %s"%ordering)
578 gross 2470 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 jfenwick 3360 :rtype: ordering in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
584 gross 2470 """
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 jfenwick 2625 :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 gross 2470 """
593     if restart == None:
594     self.__restart=restart
595     else:
596     restart=int(restart)
597     if restart<1:
598 caltinay 3789 raise ValueError("restart must be positive.")
599 gross 2470 self.__restart=restart
600 gross 2474
601 gross 2470 def getRestart(self):
602     """
603     Returns the number of iterations steps after which GMRES is performing a restart.
604 jfenwick 2625 If ``None`` is returned no restart is performed.
605 gross 2470
606 jfenwick 2625 :rtype: ``int`` or ``None``
607 gross 2470 """
608     if self.__restart < 0:
609     return None
610     else:
611     return self.__restart
612 gross 2474 def _getRestartForC(self):
613     r=self.getRestart()
614     if r == None:
615     return -1
616     else:
617     return r
618 gross 3402
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 caltinay 3789 raise ValueError("Diagonal dominance threshold must be between 0 and 1.")
629 gross 3402 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 gross 2470 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 jfenwick 2625 :param truncation: truncation
645     :type truncation: ``int``
646 gross 2470 """
647     truncation=int(truncation)
648     if truncation<1:
649 caltinay 3789 raise ValueError("truncation must be positive.")
650 gross 2470 self.__truncation=truncation
651     def getTruncation(self):
652     """
653     Returns the number of residuals in GMRES to be stored for orthogonalization
654    
655 jfenwick 2625 :rtype: ``int``
656 gross 2470 """
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 jfenwick 2625 :param iter_max: maximum number of inner iterations
663     :type iter_max: ``int``
664 gross 2470 """
665     iter_max=int(iter_max)
666     if iter_max<1:
667 caltinay 3789 raise ValueError("maximum number of inner iteration must be positive.")
668 gross 2470 self.__inner_iter_max=iter_max
669     def getInnerIterMax(self):
670     """
671     Returns maximum number of inner iteration steps
672    
673 jfenwick 2625 :rtype: ``int``
674 gross 2470 """
675     return self.__inner_iter_max
676 gross 2479 def setIterMax(self,iter_max=100000):
677 gross 2470 """
678     Sets the maximum number of iteration steps
679    
680 jfenwick 2625 :param iter_max: maximum number of iteration steps
681     :type iter_max: ``int``
682 gross 2470 """
683     iter_max=int(iter_max)
684     if iter_max<1:
685 caltinay 3789 raise ValueError("maximum number of iteration steps must be positive.")
686 gross 2470 self.__iter_max=iter_max
687     def getIterMax(self):
688     """
689     Returns maximum number of iteration steps
690    
691 jfenwick 2625 :rtype: ``int``
692 gross 2470 """
693     return self.__iter_max
694 gross 3283 def setLevelMax(self,level_max=100):
695 gross 2470 """
696     Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
697    
698 jfenwick 2625 :param level_max: maximum number of levels
699     :type level_max: ``int``
700 gross 2470 """
701     level_max=int(level_max)
702     if level_max<0:
703 caltinay 3789 raise ValueError("maximum number of coarsening levels must be non-negative.")
704 gross 2470 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 jfenwick 2625 :rtype: ``int``
710 gross 2470 """
711     return self.__level_max
712 caltinay 3517 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 artak 2688 def setCoarseningThreshold(self,theta=0.25):
729 gross 2470 """
730     Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
731    
732 jfenwick 2625 :param theta: threshold for coarsening
733     :type theta: positive ``float``
734 gross 2470 """
735     theta=float(theta)
736     if theta<0 or theta>1:
737 caltinay 3789 raise ValueError("threshold must be non-negative and less or equal 1.")
738 gross 2470 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 jfenwick 2625 :rtype: ``float``
744 gross 2470 """
745     return self.__coarsening_threshold
746 gross 3435 def setNumSweeps(self,sweeps=1):
747 gross 2470 """
748     Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
749    
750 jfenwick 2625 :param sweeps: number of sweeps
751     :type sweeps: positive ``int``
752 gross 2470 """
753     sweeps=int(sweeps)
754     if sweeps<1:
755 caltinay 3789 raise ValueError("number of sweeps must be positive.")
756 gross 2470 self.__sweeps=sweeps
757     def getNumSweeps(self):
758     """
759     Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
760    
761 jfenwick 2625 :rtype: ``int``
762 gross 2470 """
763     return self.__sweeps
764 plaub 3436 def setNumPreSweeps(self,sweeps=1):
765 gross 2470 """
766     Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
767    
768 jfenwick 2625 :param sweeps: number of sweeps
769     :type sweeps: positive ``int``
770 gross 2470 """
771     sweeps=int(sweeps)
772     if sweeps<1:
773 caltinay 3789 raise ValueError("number of sweeps must be positive.")
774 gross 2470 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 jfenwick 2625 :rtype: ``int``
780 gross 2470 """
781     return self.__pre_sweeps
782 plaub 3436 def setNumPostSweeps(self,sweeps=1):
783 gross 2470 """
784     Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
785    
786 jfenwick 2625 :param sweeps: number of sweeps
787     :type sweeps: positive ``int``
788 gross 2470 """
789     sweeps=int(sweeps)
790     if sweeps<1:
791 caltinay 3789 raise ValueError("number of sweeps must be positive.")
792 gross 2470 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 jfenwick 2625 :rtype: ``int``
798 gross 2470 """
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 jfenwick 2625 :param rtol: relative tolerance
806     :type rtol: non-negative ``float``
807 gross 2470 """
808     rtol=float(rtol)
809     if rtol<0 or rtol>1:
810 caltinay 3789 raise ValueError("tolerance must be non-negative and less or equal 1.")
811 gross 2470 self.__tolerance=rtol
812     def getTolerance(self):
813     """
814     Returns the relative tolerance for the solver
815    
816 jfenwick 2625 :rtype: ``float``
817 gross 2470 """
818     return self.__tolerance
819     def setAbsoluteTolerance(self,atol=0.):
820     """
821     Sets the absolute tolerance for the solver
822    
823 jfenwick 2625 :param atol: absolute tolerance
824     :type atol: non-negative ``float``
825 gross 2470 """
826     atol=float(atol)
827     if atol<0:
828 caltinay 3789 raise ValueError("tolerance must be non-negative.")
829 gross 2470 self.__absolute_tolerance=atol
830     def getAbsoluteTolerance(self):
831     """
832     Returns the absolute tolerance for the solver
833    
834 jfenwick 2625 :rtype: ``float``
835 gross 2470 """
836     return self.__absolute_tolerance
837    
838     def setInnerTolerance(self,rtol=0.9):
839     """
840 jfenwick 2625 Sets the relative tolerance for an inner iteration scheme for instance on the coarsest level in a multi-level scheme.
841 gross 2470
842 jfenwick 2625 :param rtol: inner relative tolerance
843     :type rtol: positive ``float``
844 gross 2470 """
845     rtol=float(rtol)
846     if rtol<=0 or rtol>1:
847 caltinay 3789 raise ValueError("tolerance must be positive and less or equal 1.")
848 gross 2470 self.__inner_tolerance=rtol
849     def getInnerTolerance(self):
850     """
851     Returns the relative tolerance for an inner iteration scheme
852    
853 jfenwick 2625 :rtype: ``float``
854 gross 2470 """
855     return self.__inner_tolerance
856     def setDropTolerance(self,drop_tol=0.01):
857     """
858     Sets the relative drop tolerance in ILUT
859    
860 jfenwick 2625 :param drop_tol: drop tolerance
861     :type drop_tol: positive ``float``
862 gross 2470 """
863     drop_tol=float(drop_tol)
864     if drop_tol<=0 or drop_tol>1:
865 caltinay 3789 raise ValueError("drop tolerance must be positive and less or equal 1.")
866 gross 2470 self.__drop_tolerance=drop_tol
867     def getDropTolerance(self):
868     """
869     Returns the relative drop tolerance in ILUT
870    
871 jfenwick 2625 :rtype: ``float``
872 gross 2470 """
873     return self.__drop_tolerance
874     def setDropStorage(self,storage=2.):
875     """
876 jfenwick 2625 Sets the maximum allowed increase in storage for ILUT. ``storage`` =2 would mean that
877 gross 2470 a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
878    
879 jfenwick 2625 :param storage: allowed storage increase
880     :type storage: ``float``
881 gross 2470 """
882     storage=float(storage)
883     if storage<1:
884 caltinay 3789 raise ValueError("allowed storage increase must be greater or equal to 1.")
885 gross 2470 self.__drop_storage=storage
886     def getDropStorage(self):
887    
888     """
889     Returns the maximum allowed increase in storage for ILUT
890    
891 jfenwick 2625 :rtype: ``float``
892 gross 2470 """
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 jfenwick 2625 :param factor: relaxation factor
899     :type factor: ``float``
900     :note: RILU with a relaxation factor 0 is identical to ILU0
901 gross 2470 """
902     factor=float(factor)
903     if factor<0:
904 caltinay 3789 raise ValueError("relaxation factor must be non-negative.")
905 gross 2470 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 jfenwick 2625 :rtype: ``float``
912 gross 2470 """
913     return self.__relaxation
914     def isSymmetric(self):
915     """
916     Checks if symmetry of the coefficient matrix is indicated.
917    
918 jfenwick 2625 :return: True if a symmetric PDE is indicated, False otherwise
919     :rtype: ``bool``
920 gross 2470 """
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 jfenwick 2625 Sets the symmetry flag for the coefficient matrix to ``flag``.
935 gross 2470
936 jfenwick 2625 :param flag: If True, the symmetry flag is set otherwise reset.
937     :type flag: ``bool``
938 gross 2470 """
939     if flag:
940     self.setSymmetryOn()
941     else:
942     self.setSymmetryOff()
943     def isVerbose(self):
944     """
945 jfenwick 2625 Returns ``True`` if the solver is expected to be verbose.
946 gross 2470
947 jfenwick 2625 :return: True if verbosity of switched on.
948     :rtype: ``bool``
949 gross 2470 """
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 gross 2474 def setVerbosity(self,verbose=False):
963 gross 2470 """
964 jfenwick 2625 Sets the verbosity flag for the solver to ``flag``.
965 gross 2470
966 jfenwick 2625 :param verbose: If ``True``, the verbosity of the solver is switched on.
967     :type verbose: ``bool``
968 gross 2470 """
969     if verbose:
970     self.setVerbosityOn()
971     else:
972     self.setVerbosityOff()
973 gross 2474
974 gross 2470 def adaptInnerTolerance(self):
975     """
976 jfenwick 2625 Returns ``True`` if the tolerance of the inner solver is selected automatically.
977     Otherwise the inner tolerance set by `setInnerTolerance` is used.
978 gross 2470
979 jfenwick 2625 :return: ``True`` if inner tolerance adaption is chosen.
980     :rtype: ``bool``
981 gross 2470 """
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 gross 3352 Sets the flag to indicate automatic selection of the inner tolerance.
997 gross 2470
998 jfenwick 2625 :param adapt: If ``True``, the inner tolerance is selected automatically.
999     :type adapt: ``bool``
1000 gross 2470 """
1001     if adapt:
1002     self.setInnerToleranceAdaptionOn()
1003     else:
1004     self.setInnerToleranceAdaptionOff()
1005    
1006     def acceptConvergenceFailure(self):
1007     """
1008 jfenwick 2625 Returns ``True`` if a failure to meet the stopping criteria within the
1009 gross 2470 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 jfenwick 2625 stopping criteria. One can use the `hasConverged` method to check if the
1013 gross 2470 last call to the solver was successful.
1014    
1015 jfenwick 2625 :return: ``True`` if a failure to achieve convergence is accepted.
1016     :rtype: ``bool``
1017 gross 2470 """
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 gross 3352 Sets the flag to indicate the acceptance of a failure of convergence.
1033 gross 2470
1034 jfenwick 2625 :param accept: If ``True``, any failure to achieve convergence is accepted.
1035     :type accept: ``bool``
1036 gross 2470 """
1037     if accept:
1038     self.setAcceptanceConvergenceFailureOn()
1039     else:
1040     self.setAcceptanceConvergenceFailureOff()
1041    
1042 gross 3094 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 gross 3352 Sets the flag to use local preconditioning to on
1056 gross 3094 """
1057     self.__use_local_preconditioner=True
1058     def setLocalPreconditionerOff(self):
1059     """
1060 gross 3352 Sets the flag to use local preconditioning to off
1061 gross 3094 """
1062     self.__use_local_preconditioner=False
1063 gross 3352
1064 gross 3094 def setLocalPreconditioner(self,use=False):
1065     """
1066 gross 3352 Sets the flag to use local preconditioning
1067 gross 3094
1068 gross 3352 :param use: If ``True``, local proconditioning on each MPI rank is applied
1069     :type use: ``bool``
1070 gross 3094 """
1071     if use:
1072     self.setUseLocalPreconditionerOn()
1073     else:
1074     self.setUseLocalPreconditionerOff()
1075 gross 3283
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 caltinay 3789 raise ValueError("sparsity must be non-negative.")
1088 gross 3283 if sparsity>1:
1089 caltinay 3789 raise ValueError("sparsity must be less than 1.")
1090 gross 3283 self.__min_sparsity=sparsity
1091 gross 3094
1092 gross 3283 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 caltinay 3789 raise ValueError("number of refinements must be non-negative.")
1113 gross 3283 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 caltinay 3789 raise ValueError("number of refinements must be non-negative.")
1133 gross 3283 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 gross 3352
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 jfenwick 3360 :param use: If ``True``,a panel is used to find unknowns in AMG coarsening
1170 gross 3352 :type use: ``bool``
1171     """
1172     if use:
1173     self.setUsePanelOn()
1174     else:
1175     self.setUsePanelOff()
1176    
1177 gross 3440 def setAMGInterpolation(self, method=None):
1178 gross 3352 """
1179 gross 3440 Set the interpolation method for the AMG preconditioner.
1180 gross 3352
1181 gross 3440 :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 gross 3352 """
1184 gross 3440 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 caltinay 3789 raise ValueError("unknown AMG interpolation method %s"%method)
1187 gross 3440 self.__amg_interpolation_method=method
1188 gross 3352
1189 gross 3440 def getAMGInterpolation(self):
1190 gross 3352 """
1191 gross 3440 Returns key of the interpolation method for the AMG preconditioner
1192 gross 3352
1193 gross 3440 :rtype: in the list `SolverOptions.CLASSIC_INTERPOLATION_WITH_FF_COUPLING`, `SolverOptions.CLASSIC_INTERPOLATION', `SolverOptions.DIRECT_INTERPOLATION`
1194 gross 3352 """
1195 gross 3440 return self.__amg_interpolation_method
1196 gross 3352
1197    
1198 jgs 148 class IllegalCoefficient(ValueError):
1199 jgs 102 """
1200 caltinay 2169 Exception that is raised if an illegal coefficient of the general or
1201     particular PDE is requested.
1202 jgs 148 """
1203 gross 1072 pass
1204 jgs 102
1205 jgs 148 class IllegalCoefficientValue(ValueError):
1206 jgs 102 """
1207 caltinay 2169 Exception that is raised if an incorrect value for a coefficient is used.
1208 jgs 148 """
1209 gross 1072 pass
1210 jgs 122
1211 gross 1072 class IllegalCoefficientFunctionSpace(ValueError):
1212     """
1213 caltinay 2169 Exception that is raised if an incorrect function space for a coefficient
1214     is used.
1215 gross 1072 """
1216    
1217 jgs 148 class UndefinedPDEError(ValueError):
1218     """
1219 caltinay 2169 Exception that is raised if a PDE is not fully defined yet.
1220 jgs 148 """
1221 gross 1072 pass
1222 jgs 102
1223 gross 1841 class PDECoef(object):
1224 jgs 102 """
1225 caltinay 2169 A class for describing a PDE coefficient.
1226 jgs 149
1227 jfenwick 2625 :cvar INTERIOR: indicator that coefficient is defined on the interior of
1228 caltinay 2169 the PDE domain
1229 jfenwick 2625 :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
1230 caltinay 2169 the PDE domain
1231 jfenwick 2625 :cvar CONTACT: indicator that coefficient is defined on the contact region
1232 caltinay 2169 within the PDE domain
1233 jfenwick 2625 :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
1234 caltinay 2169 interior of the PDE domain using a reduced
1235     integration order
1236 jfenwick 2625 :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
1237 caltinay 2169 boundary of the PDE domain using a reduced
1238     integration order
1239 jfenwick 2625 :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
1240 caltinay 2169 region within the PDE domain using a reduced
1241     integration order
1242 jfenwick 2625 :cvar SOLUTION: indicator that coefficient is defined through a solution of
1243 caltinay 2169 the PDE
1244 jfenwick 2625 :cvar REDUCED: indicator that coefficient is defined through a reduced
1245 caltinay 2169 solution of the PDE
1246 caltinay 3609 :cvar DIRACDELTA: indicator that coefficient is defined as Dirac delta functions
1247 jfenwick 2625 :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
1248 caltinay 2169 defined by the number of PDE equations
1249 jfenwick 2625 :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
1250 caltinay 2169 defined by the number of PDE solutions
1251 jfenwick 2625 :cvar BY_DIM: indicator that the dimension of the coefficient shape is
1252 caltinay 2169 defined by the spatial dimension
1253 jfenwick 2625 :cvar OPERATOR: indicator that the the coefficient alters the operator of
1254 caltinay 2169 the PDE
1255 jfenwick 2625 :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
1256 caltinay 2169 hand side of the PDE
1257 jfenwick 2625 :cvar BOTH: indicator that the the coefficient alters the operator as well
1258 caltinay 2169 as the right hand side of the PDE
1259 jgs 149
1260 jgs 102 """
1261     INTERIOR=0
1262     BOUNDARY=1
1263     CONTACT=2
1264 jgs 149 SOLUTION=3
1265 jgs 150 REDUCED=4
1266 jgs 149 BY_EQUATION=5
1267     BY_SOLUTION=6
1268     BY_DIM=7
1269     OPERATOR=10
1270     RIGHTHANDSIDE=11
1271     BOTH=12
1272 gross 1072 INTERIOR_REDUCED=13
1273     BOUNDARY_REDUCED=14
1274     CONTACT_REDUCED=15
1275 caltinay 3609 DIRACDELTA=16
1276 jgs 149
1277 gross 1072 def __init__(self, where, pattern, altering):
1278 jgs 102 """
1279 caltinay 2169 Initialises a PDE coefficient type.
1280 jgs 151
1281 jfenwick 2625 :param where: describes where the coefficient lives
1282     :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1283     `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1284 caltinay 3609 `CONTACT_REDUCED`, 'DIRACDELTA'
1285 jfenwick 2625 :param pattern: describes the shape of the coefficient and how the shape
1286 caltinay 2169 is built for a given spatial dimension and numbers of
1287     equations and solutions in then PDE. For instance,
1288 jfenwick 2625 (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1289 caltinay 2169 rank 3 coefficient which is instantiated as shape (3,2,2)
1290     in case of three equations and two solution components
1291     on a 2-dimensional domain. In the case of single equation
1292     and a single solution component the shape components
1293 jfenwick 2625 marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1294 caltinay 2169 In this case the example would be read as (2,).
1295 jfenwick 2625 :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1296     :param altering: indicates what part of the PDE is altered if the
1297 caltinay 2169 coefficient is altered
1298 jfenwick 2625 :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1299 jgs 102 """
1300 gross 1841 super(PDECoef, self).__init__()
1301 jgs 102 self.what=where
1302     self.pattern=pattern
1303     self.altering=altering
1304 jgs 108 self.resetValue()
1305 jgs 102
1306 jgs 108 def resetValue(self):
1307     """
1308 caltinay 2169 Resets the coefficient value to the default.
1309 jgs 108 """
1310     self.value=escript.Data()
1311    
1312 jgs 149 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1313 jgs 102 """
1314 jfenwick 2625 Returns the `FunctionSpace` of the coefficient.
1315 jgs 102
1316 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1317     :type domain: `Domain`
1318     :param reducedEquationOrder: True to indicate that reduced order is used
1319 caltinay 2169 to represent the equation
1320 jfenwick 2625 :type reducedEquationOrder: ``bool``
1321     :param reducedSolutionOrder: True to indicate that reduced order is used
1322 caltinay 2169 to represent the solution
1323 jfenwick 2625 :type reducedSolutionOrder: ``bool``
1324     :return: `FunctionSpace` of the coefficient
1325     :rtype: `FunctionSpace`
1326 jgs 102 """
1327 jgs 151 if self.what==self.INTERIOR:
1328 jgs 149 return escript.Function(domain)
1329 gross 1072 elif self.what==self.INTERIOR_REDUCED:
1330     return escript.ReducedFunction(domain)
1331 jgs 151 elif self.what==self.BOUNDARY:
1332 jgs 149 return escript.FunctionOnBoundary(domain)
1333 gross 1072 elif self.what==self.BOUNDARY_REDUCED:
1334     return escript.ReducedFunctionOnBoundary(domain)
1335 jgs 151 elif self.what==self.CONTACT:
1336 jgs 149 return escript.FunctionOnContactZero(domain)
1337 gross 1072 elif self.what==self.CONTACT_REDUCED:
1338     return escript.ReducedFunctionOnContactZero(domain)
1339 caltinay 3609 elif self.what==self.DIRACDELTA:
1340     return escript.DiracDeltaFunctions(domain)
1341 jgs 151 elif self.what==self.SOLUTION:
1342 jgs 149 if reducedEquationOrder and reducedSolutionOrder:
1343     return escript.ReducedSolution(domain)
1344     else:
1345     return escript.Solution(domain)
1346 jgs 151 elif self.what==self.REDUCED:
1347     return escript.ReducedSolution(domain)
1348 jgs 102
1349 jgs 108 def getValue(self):
1350     """
1351 caltinay 2169 Returns the value of the coefficient.
1352 jgs 149
1353 jfenwick 2625 :return: value of the coefficient
1354     :rtype: `Data`
1355 jgs 108 """
1356     return self.value
1357 jgs 148
1358 jgs 149 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
1359 jgs 108 """
1360 caltinay 2169 Sets the value of the coefficient to a new value.
1361 jgs 149
1362 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1363     :type domain: `Domain`
1364     :param numEquations: number of equations of the PDE
1365     :type numEquations: ``int``
1366     :param numSolutions: number of components of the PDE solution
1367     :type numSolutions: ``int``
1368     :param reducedEquationOrder: True to indicate that reduced order is used
1369 caltinay 2169 to represent the equation
1370 jfenwick 2625 :type reducedEquationOrder: ``bool``
1371     :param reducedSolutionOrder: True to indicate that reduced order is used
1372 caltinay 2169 to represent the solution
1373 jfenwick 2625 :type reducedSolutionOrder: ``bool``
1374     :param newValue: number of components of the PDE solution
1375     :type newValue: any object that can be converted into a
1376     `Data` object with the appropriate shape
1377     and `FunctionSpace`
1378     :raise IllegalCoefficientValue: if the shape of the assigned value does
1379 caltinay 2169 not match the shape of the coefficient
1380 jfenwick 2625 :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1381 caltinay 2169 to appropriate function space
1382 jgs 108 """
1383 jgs 148 if newValue==None:
1384     newValue=escript.Data()
1385     elif isinstance(newValue,escript.Data):
1386     if not newValue.isEmpty():
1387 gross 1072 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
1388     try:
1389     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1390     except:
1391 caltinay 3789 raise IllegalCoefficientFunctionSpace("Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain))
1392 jgs 148 else:
1393 jgs 149 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1394 jgs 148 if not newValue.isEmpty():
1395     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1396 caltinay 3789 raise IllegalCoefficientValue("Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape()))
1397 jgs 108 self.value=newValue
1398 jgs 148
1399 jgs 102 def isAlteringOperator(self):
1400     """
1401 caltinay 2169 Checks if the coefficient alters the operator of the PDE.
1402 jgs 149
1403 jfenwick 2625 :return: True if the operator of the PDE is changed when the
1404 caltinay 2169 coefficient is changed
1405 jfenwick 2625 :rtype: ``bool``
1406 caltinay 2169 """
1407 jgs 102 if self.altering==self.OPERATOR or self.altering==self.BOTH:
1408     return not None
1409     else:
1410     return None
1411    
1412     def isAlteringRightHandSide(self):
1413     """
1414 caltinay 2169 Checks if the coefficient alters the right hand side of the PDE.
1415 jgs 149
1416 jfenwick 2625 :rtype: ``bool``
1417     :return: True if the right hand side of the PDE is changed when the
1418     coefficient is changed, ``None`` otherwise.
1419 caltinay 2169 """
1420 jgs 102 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1421     return not None
1422     else:
1423     return None
1424    
1425 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
1426 jgs 102 """
1427 caltinay 2169 Tries to estimate the number of equations and number of solutions if
1428     the coefficient has the given shape.
1429 jgs 102
1430 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1431     :type domain: `Domain`
1432     :param shape: suggested shape of the coefficient
1433     :type shape: ``tuple`` of ``int`` values
1434     :return: the number of equations and number of solutions of the PDE if
1435 caltinay 2169 the coefficient has given shape. If no appropriate numbers
1436 jfenwick 2625 could be identified, ``None`` is returned
1437     :rtype: ``tuple`` of two ``int`` values or ``None``
1438 jgs 102 """
1439 jgs 148 dim=domain.getDim()
1440 jgs 102 if len(shape)>0:
1441     num=max(shape)+1
1442     else:
1443     num=1
1444     search=[]
1445 jgs 149 if self.definesNumEquation() and self.definesNumSolutions():
1446     for u in range(num):
1447     for e in range(num):
1448     search.append((e,u))
1449     search.sort(self.__CompTuple2)
1450     for item in search:
1451 jgs 148 s=self.getShape(domain,item[0],item[1])
1452 jgs 102 if len(s)==0 and len(shape)==0:
1453     return (1,1)
1454     else:
1455     if s==shape: return item
1456 jgs 149 elif self.definesNumEquation():
1457     for e in range(num,0,-1):
1458     s=self.getShape(domain,e,0)
1459     if len(s)==0 and len(shape)==0:
1460     return (1,None)
1461     else:
1462     if s==shape: return (e,None)
1463    
1464     elif self.definesNumSolutions():
1465     for u in range(num,0,-1):
1466     s=self.getShape(domain,0,u)
1467     if len(s)==0 and len(shape)==0:
1468     return (None,1)
1469     else:
1470     if s==shape: return (None,u)
1471 jgs 102 return None
1472 caltinay 2169
1473 jgs 149 def definesNumSolutions(self):
1474     """
1475 caltinay 2169 Checks if the coefficient allows to estimate the number of solution
1476     components.
1477 jgs 102
1478 jfenwick 2625 :return: True if the coefficient allows an estimate of the number of
1479 caltinay 2169 solution components, False otherwise
1480 jfenwick 2625 :rtype: ``bool``
1481 jgs 149 """
1482     for i in self.pattern:
1483     if i==self.BY_SOLUTION: return True
1484     return False
1485    
1486     def definesNumEquation(self):
1487     """
1488 caltinay 2169 Checks if the coefficient allows to estimate the number of equations.
1489 jgs 149
1490 jfenwick 2625 :return: True if the coefficient allows an estimate of the number of
1491 caltinay 2169 equations, False otherwise
1492 jfenwick 2625 :rtype: ``bool``
1493 jgs 149 """
1494     for i in self.pattern:
1495     if i==self.BY_EQUATION: return True
1496     return False
1497    
1498     def __CompTuple2(self,t1,t2):
1499     """
1500 caltinay 2169 Compares two tuples of possible number of equations and number of
1501     solutions.
1502 jgs 149
1503 jfenwick 2625 :param t1: the first tuple
1504     :param t2: the second tuple
1505     :return: 0, 1, or -1
1506 jgs 149 """
1507    
1508     dif=t1[0]+t1[1]-(t2[0]+t2[1])
1509     if dif<0: return 1
1510     elif dif>0: return -1
1511     else: return 0
1512    
1513 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
1514 jgs 149 """
1515 caltinay 2169 Builds the required shape of the coefficient.
1516 jgs 102
1517 jfenwick 2625 :param domain: domain on which the PDE uses the coefficient
1518     :type domain: `Domain`
1519     :param numEquations: number of equations of the PDE
1520     :type numEquations: ``int``
1521     :param numSolutions: number of components of the PDE solution
1522     :type numSolutions: ``int``
1523     :return: shape of the coefficient
1524     :rtype: ``tuple`` of ``int`` values
1525 jgs 149 """
1526     dim=domain.getDim()
1527     s=()
1528     for i in self.pattern:
1529     if i==self.BY_EQUATION:
1530 jgs 148 if numEquations>1: s=s+(numEquations,)
1531 jgs 149 elif i==self.BY_SOLUTION:
1532 jgs 148 if numSolutions>1: s=s+(numSolutions,)
1533 jgs 102 else:
1534     s=s+(dim,)
1535 jgs 149 return s
1536 jgs 102
1537 gross 2470 #====================================================================================================================
1538    
1539 gross 1841 class LinearProblem(object):
1540 jgs 102 """
1541 caltinay 2169 This is the base class to define a general linear PDE-type problem for
1542 jfenwick 2625 for an unknown function *u* on a given domain defined through a
1543     `Domain` object. The problem can be given as a single
1544 caltinay 2169 equation or as a system of equations.
1545 jgs 102
1546 caltinay 2169 The class assumes that some sort of assembling process is required to form
1547     a problem of the form
1548 jgs 151
1549 jfenwick 2625 *L u=f*
1550 jgs 102
1551 jfenwick 2625 where *L* is an operator and *f* is the right hand side. This operator
1552     problem will be solved to get the unknown *u*.
1553 caltinay 2169
1554 jgs 102 """
1555 jgs 148 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1556 jgs 102 """
1557 caltinay 2169 Initializes a linear problem.
1558 jgs 102
1559 jfenwick 2625 :param domain: domain of the PDE
1560     :type domain: `Domain`
1561     :param numEquations: number of equations. If ``None`` the number of
1562 caltinay 2169 equations is extracted from the coefficients.
1563 jfenwick 2625 :param numSolutions: number of solution components. If ``None`` the number
1564 caltinay 2169 of solution components is extracted from the
1565     coefficients.
1566 jfenwick 2625 :param debug: if True debug information is printed
1567 jgs 148
1568 jgs 102 """
1569 gross 1841 super(LinearProblem, self).__init__()
1570 jgs 102
1571 jgs 148 self.__debug=debug
1572 jgs 104 self.__domain=domain
1573     self.__numEquations=numEquations
1574     self.__numSolutions=numSolutions
1575 gross 1841 self.__altered_coefficients=False
1576 jgs 149 self.__reduce_equation_order=False
1577     self.__reduce_solution_order=False
1578 gross 2474 self.__sym=False
1579     self.setSolverOptions()
1580 gross 1841 self.__is_RHS_valid=False
1581     self.__is_operator_valid=False
1582     self.__COEFFICIENTS={}
1583 gross 2474 self.__solution_rtol=1.e99
1584     self.__solution_atol=1.e99
1585     self.setSymmetryOff()
1586 gross 1841 # initialize things:
1587     self.resetAllCoefficients()
1588 gross 2474 self.initializeSystem()
1589 caltinay 2169 # ==========================================================================
1590 jgs 148 # general stuff:
1591 caltinay 2169 # ==========================================================================
1592 jgs 148 def __str__(self):
1593 jgs 149 """
1594 caltinay 2169 Returns a string representation of the PDE.
1595 jgs 149
1596 jfenwick 2625 :return: a simple representation of the PDE
1597     :rtype: ``str``
1598 jgs 149 """
1599 gross 1841 return "<LinearProblem %d>"%id(self)
1600 caltinay 2169 # ==========================================================================
1601 jgs 148 # debug :
1602 caltinay 2169 # ==========================================================================
1603 jgs 148 def setDebugOn(self):
1604 jgs 149 """
1605 caltinay 2169 Switches debug output on.
1606 jgs 108 """
1607 jgs 148 self.__debug=not None
1608    
1609     def setDebugOff(self):
1610 jgs 108 """
1611 caltinay 2169 Switches debug output off.
1612 jgs 148 """
1613     self.__debug=None
1614 jgs 108
1615 gross 2100 def setDebug(self, flag):
1616     """
1617 jfenwick 2625 Switches debug output on if ``flag`` is True otherwise it is switched off.
1618 caltinay 2169
1619 jfenwick 2625 :param flag: desired debug status
1620     :type flag: ``bool``
1621 gross 2100 """
1622     if flag:
1623     self.setDebugOn()
1624     else:
1625     self.setDebugOff()
1626    
1627 jgs 148 def trace(self,text):
1628     """
1629 caltinay 2169 Prints the text message if debug mode is switched on.
1630    
1631 jfenwick 2625 :param text: message to be printed
1632     :type text: ``string``
1633 jgs 102 """
1634 caltinay 3789 if self.__debug: print(("%s: %s"%(str(self),text)))
1635 jgs 102
1636 caltinay 2169 # ==========================================================================
1637 jgs 148 # some service functions:
1638 caltinay 2169 # ==========================================================================
1639 gross 1841 def introduceCoefficients(self,**coeff):
1640     """
1641 caltinay 2169 Introduces new coefficients into the problem.
1642 gross 1841
1643 jfenwick 2625 Use:
1644 gross 1841
1645 caltinay 2169 p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1646 gross 1841
1647 jfenwick 2625 to introduce the coefficients *A* and *B*.
1648 gross 1841 """
1649 caltinay 3789 for name, type in list(coeff.items()):
1650 gross 1841 if not isinstance(type,PDECoef):
1651 caltinay 3789 raise ValueError("coefficient %s has no type."%name)
1652 gross 1841 self.__COEFFICIENTS[name]=type
1653     self.__COEFFICIENTS[name].resetValue()
1654     self.trace("coefficient %s has been introduced."%name)
1655    
1656 jgs 148 def getDomain(self):
1657 jgs 102 """
1658 caltinay 2169 Returns the domain of the PDE.
1659 jgs 102
1660 jfenwick 2625 :return: the domain of the PDE
1661     :rtype: `Domain`
1662 jgs 108 """
1663 jgs 148 return self.__domain
1664 gross 2535 def getDomainStatus(self):
1665     """
1666     Return the status indicator of the domain
1667     """
1668     return self.getDomain().getStatus()
1669 jgs 122
1670 gross 2535 def getSystemStatus(self):
1671     """
1672     Return the domain status used to build the current system
1673     """
1674     return self.__system_status
1675     def setSystemStatus(self,status=None):
1676     """
1677 jfenwick 2625 Sets the system status to ``status`` if ``status`` is not present the
1678 gross 2535 current status of the domain is used.
1679     """
1680     if status == None:
1681     self.__system_status=self.getDomainStatus()
1682     else:
1683     self.__system_status=status
1684    
1685 jgs 148 def getDim(self):
1686 jgs 108 """
1687 caltinay 2169 Returns the spatial dimension of the PDE.
1688 jgs 108
1689 jfenwick 2625 :return: the spatial dimension of the PDE domain
1690     :rtype: ``int``
1691 jgs 148 """
1692     return self.getDomain().getDim()
1693 jgs 102
1694 jgs 148 def getNumEquations(self):
1695     """
1696 caltinay 2169 Returns the number of equations.
1697 jgs 102
1698 jfenwick 2625 :return: the number of equations
1699     :rtype: ``int``
1700     :raise UndefinedPDEError: if the number of equations is not specified yet
1701 jgs 148 """
1702     if self.__numEquations==None:
1703 gross 1859 if self.__numSolutions==None:
1704 caltinay 3789 raise UndefinedPDEError("Number of equations is undefined. Please specify argument numEquations.")
1705 gross 1859 else:
1706     self.__numEquations=self.__numSolutions
1707     return self.__numEquations
1708 jgs 147
1709 jgs 148 def getNumSolutions(self):
1710     """
1711 caltinay 2169 Returns the number of unknowns.
1712 jgs 147
1713 jfenwick 2625 :return: the number of unknowns
1714     :rtype: ``int``
1715     :raise UndefinedPDEError: if the number of unknowns is not specified yet
1716 jgs 148 """
1717     if self.__numSolutions==None:
1718 gross 1859 if self.__numEquations==None:
1719 caltinay 3789 raise UndefinedPDEError("Number of solution is undefined. Please specify argument numSolutions.")
1720 gross 1859 else:
1721     self.__numSolutions=self.__numEquations
1722     return self.__numSolutions
1723 jgs 148
1724 jgs 149 def reduceEquationOrder(self):
1725     """
1726 caltinay 2169 Returns the status of order reduction for the equation.
1727 jgs 149
1728 jfenwick 2625 :return: True if reduced interpolation order is used for the
1729 caltinay 2169 representation of the equation, False otherwise
1730 jfenwick 2625 :rtype: `bool`
1731 jgs 149 """
1732     return self.__reduce_equation_order
1733    
1734     def reduceSolutionOrder(self):
1735     """
1736 caltinay 2169 Returns the status of order reduction for the solution.
1737 jgs 149
1738 jfenwick 2625 :return: True if reduced interpolation order is used for the
1739 caltinay 2169 representation of the solution, False otherwise
1740 jfenwick 2625 :rtype: `bool`
1741 jgs 149 """
1742     return self.__reduce_solution_order
1743 jgs 151
1744 jgs 108 def getFunctionSpaceForEquation(self):
1745     """
1746 jfenwick 2625 Returns the `FunctionSpace` used to discretize
1747 caltinay 2169 the equation.
1748 jgs 148
1749 jfenwick 2625 :return: representation space of equation
1750     :rtype: `FunctionSpace`
1751 jgs 108 """
1752 jgs 149 if self.reduceEquationOrder():
1753     return escript.ReducedSolution(self.getDomain())
1754     else:
1755     return escript.Solution(self.getDomain())
1756 jgs 108
1757     def getFunctionSpaceForSolution(self):
1758     """
1759 jfenwick 2625 Returns the `FunctionSpace` used to represent
1760 caltinay 2169 the solution.
1761 jgs 148
1762 jfenwick 2625 :return: representation space of solution
1763     :rtype: `FunctionSpace`
1764 jgs 108 """
1765 jgs 149 if self.reduceSolutionOrder():
1766     return escript.ReducedSolution(self.getDomain())
1767     else:
1768     return escript.Solution(self.getDomain())
1769 jgs 108
1770 caltinay 2169 # ==========================================================================
1771 jgs 148 # solver settings:
1772 caltinay 2169 # ==========================================================================
1773 gross 2474 def setSolverOptions(self,options=None):
1774 jgs 102 """
1775 gross 2474 Sets the solver options.
1776 jgs 148
1777 jfenwick 2625 :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1778     :type options: `SolverOptions` or ``None``
1779     :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1780 jgs 102 """
1781 gross 2474 if options==None:
1782     self.__solver_options=SolverOptions()
1783     elif isinstance(options, SolverOptions):
1784     self.__solver_options=options
1785     else:
1786 caltinay 3789 raise ValueError("options must be a SolverOptions object.")
1787 gross 2474 self.__solver_options.setSymmetry(self.__sym)
1788    
1789     def getSolverOptions(self):
1790 jgs 148 """
1791 gross 2474 Returns the solver options
1792    
1793 jfenwick 2625 :rtype: `SolverOptions`
1794 jgs 148 """
1795 gross 2474 self.__solver_options.setSymmetry(self.__sym)
1796     return self.__solver_options
1797    
1798 jgs 148 def isUsingLumping(self):
1799     """
1800 caltinay 2169 Checks if matrix lumping is the current solver method.
1801 jgs 148
1802 jfenwick 2625 :return: True if the current solver method is lumping
1803     :rtype: ``bool``
1804 jgs 148 """
1805 gross 3379 return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
1806 caltinay 2169 # ==========================================================================
1807 jgs 148 # symmetry flag:
1808 caltinay 2169 # ==========================================================================
1809     def isSymmetric(self):
1810 jgs 102 """
1811 caltinay 2169 Checks if symmetry is indicated.
1812 jgs 149
1813 jfenwick 2625 :return: True if a symmetric PDE is indicated, False otherwise
1814     :rtype: ``bool``
1815     :note: the method is equivalent to use getSolverOptions().isSymmetric()
1816 jgs 102 """
1817 gross 2474 self.getSolverOptions().isSymmetric()
1818 jgs 102
1819 caltinay 2169 def setSymmetryOn(self):
1820 jgs 102 """
1821 gross 2474 Sets the symmetry flag.
1822 jfenwick 2625 :note: The method overwrites the symmetry flag set by the solver options
1823 jgs 102 """
1824 gross 2474 self.__sym=True
1825     self.getSolverOptions().setSymmetryOn()
1826 jgs 102
1827 caltinay 2169 def setSymmetryOff(self):
1828 jgs 102 """
1829 caltinay 2169 Clears the symmetry flag.
1830 jfenwick 2625 :note: The method overwrites the symmetry flag set by the solver options
1831 jgs 102 """
1832 gross 2474 self.__sym=False
1833     self.getSolverOptions().setSymmetryOff()
1834 jgs 102
1835 gross 2474 def setSymmetry(self,flag=False):
1836 jgs 148 """
1837 jfenwick 2625 Sets the symmetry flag to ``flag``.
1838 jgs 149
1839 jfenwick 2625 :param flag: If True, the symmetry flag is set otherwise reset.
1840     :type flag: ``bool``
1841     :note: The method overwrites the symmetry flag set by the solver options
1842 jgs 148 """
1843 gross 2474 self.getSolverOptions().setSymmetry(flag)
1844 caltinay 2169 # ==========================================================================
1845 jgs 148 # function space handling for the equation as well as the solution
1846 caltinay 2169 # ==========================================================================
1847     def setReducedOrderOn(self):
1848 jgs 102 """
1849 caltinay 2169 Switches reduced order on for solution and equation representation.
1850 jgs 149
1851 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1852 caltinay 2169 been set
1853 jgs 102 """
1854     self.setReducedOrderForSolutionOn()
1855     self.setReducedOrderForEquationOn()
1856    
1857 caltinay 2169 def setReducedOrderOff(self):
1858 jgs 102 """
1859 caltinay 2169 Switches reduced order off for solution and equation representation
1860 jgs 149
1861 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1862 caltinay 2169 been set
1863 jgs 102 """
1864     self.setReducedOrderForSolutionOff()
1865     self.setReducedOrderForEquationOff()
1866    
1867 caltinay 2169 def setReducedOrderTo(self,flag=False):
1868 jgs 102 """
1869 caltinay 2169 Sets order reduction state for both solution and equation representation
1870     according to flag.
1871    
1872 jfenwick 2625 :param flag: if True, the order reduction is switched on for both solution
1873 caltinay 2169 and equation representation, otherwise or if flag is not
1874     present order reduction is switched off
1875 jfenwick 2625 :type flag: ``bool``
1876     :raise RuntimeError: if order reduction is altered after a coefficient has
1877 caltinay 2169 been set
1878 jgs 102 """
1879     self.setReducedOrderForSolutionTo(flag)
1880     self.setReducedOrderForEquationTo(flag)
1881    
1882 jgs 148
1883 caltinay 2169 def setReducedOrderForSolutionOn(self):
1884 jgs 102 """
1885 caltinay 2169 Switches reduced order on for solution representation.
1886 jgs 149
1887 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1888 caltinay 2169 been set
1889 jgs 102 """
1890 jgs 149 if not self.__reduce_solution_order:
1891     if self.__altered_coefficients:
1892 caltinay 3789 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1893 caltinay 2169 self.trace("Reduced order is used for solution representation.")
1894 jgs 149 self.__reduce_solution_order=True
1895 gross 1841 self.initializeSystem()
1896 jgs 102
1897 caltinay 2169 def setReducedOrderForSolutionOff(self):
1898 jgs 102 """
1899 caltinay 2169 Switches reduced order off for solution representation
1900 jgs 149
1901 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1902 caltinay 2169 been set.
1903 jgs 102 """
1904 jgs 149 if self.__reduce_solution_order:
1905     if self.__altered_coefficients:
1906 caltinay 3789 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1907 jgs 148 self.trace("Full order is used to interpolate solution.")
1908 jgs 149 self.__reduce_solution_order=False
1909 gross 1841 self.initializeSystem()
1910 jgs 102
1911 caltinay 2169 def setReducedOrderForSolutionTo(self,flag=False):
1912 jgs 102 """
1913 caltinay 2169 Sets order reduction state for solution representation according to flag.
1914 jgs 102
1915 jfenwick 2625 :param flag: if flag is True, the order reduction is switched on for
1916 caltinay 2169 solution representation, otherwise or if flag is not present
1917     order reduction is switched off
1918 jfenwick 2625 :type flag: ``bool``
1919     :raise RuntimeError: if order reduction is altered after a coefficient has
1920 caltinay 2169 been set
1921 jgs 102 """
1922     if flag:
1923     self.setReducedOrderForSolutionOn()
1924     else:
1925     self.setReducedOrderForSolutionOff()
1926 jgs 148
1927 caltinay 2169 def setReducedOrderForEquationOn(self):
1928 jgs 102 """
1929 caltinay 2169 Switches reduced order on for equation representation.
1930 jgs 149
1931 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1932 caltinay 2169 been set
1933 jgs 102 """
1934 jgs 149 if not self.__reduce_equation_order:
1935     if self.__altered_coefficients:
1936 caltinay 3789 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1937 jgs 148 self.trace("Reduced order is used for test functions.")
1938 jgs 149 self.__reduce_equation_order=True
1939 gross 1841 self.initializeSystem()
1940 jgs 102
1941 caltinay 2169 def setReducedOrderForEquationOff(self):
1942 jgs 102 """
1943 caltinay 2169 Switches reduced order off for equation representation.
1944 jgs 149
1945 jfenwick 2625 :raise RuntimeError: if order reduction is altered after a coefficient has
1946 caltinay 2169 been set
1947 jgs 102 """
1948 jgs 149 if self.__reduce_equation_order:
1949     if self.__altered_coefficients:
1950 caltinay 3789 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1951 jgs 148 self.trace("Full order is used for test functions.")
1952 jgs 149 self.__reduce_equation_order=False
1953 gross 1841 self.initializeSystem()
1954 jgs 102
1955 caltinay 2169 def setReducedOrderForEquationTo(self,flag=False):
1956 jgs 102 """
1957 caltinay 2169 Sets order reduction state for equation representation according to flag.
1958 jgs 102
1959 jfenwick 2625 :param flag: if flag is True, the order reduction is switched on for
1960 caltinay 2169 equation representation, otherwise or if flag is not present
1961     order reduction is switched off
1962 jfenwick 2625 :type flag: ``bool``
1963     :raise RuntimeError: if order reduction is altered after a coefficient has
1964 caltinay 2169 been set
1965 jgs 102 """
1966     if flag:
1967     self.setReducedOrderForEquationOn()
1968     else:
1969     self.setReducedOrderForEquationOff()
1970 gross 2474 def getOperatorType(self):
1971 gross 1841 """
1972 caltinay 2169 Returns the current system type.
1973 gross 1841 """
1974 gross 2474 return self.__operator_type
1975 jgs 148
1976 gross 1841 def checkSymmetricTensor(self,name,verbose=True):
1977     """
1978 caltinay 2169 Tests a coefficient for symmetry.
1979 jgs 148
1980 jfenwick 2625 :param name: name of the coefficient
1981     :type name: ``str``
1982     :param verbose: if set to True or not present a report on coefficients
1983 caltinay 2169 which break the symmetry is printed.
1984 jfenwick 2625 :type verbose: ``bool``
1985     :return: True if coefficient ``name`` is symmetric
1986     :rtype: ``bool``
1987 gross 1841 """
1988 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
1989 gross 1841 A=self.getCoefficient(name)
1990     verbose=verbose or self.__debug
1991     out=True
1992     if not A.isEmpty():
1993 gross 2474 tol=util.Lsup(A)*SMALL_TOLERANCE
1994 gross 1841 s=A.getShape()
1995     if A.getRank() == 4:
1996     if s[0]==s[2] and s[1] == s[3]:
1997     for i in range(s[0]):
1998     for j in range(s[1]):
1999     for k in range(s[2]):
2000     for l in range(s[3]):
2001     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
2002 caltinay 3789 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)))
2003 gross 1841 out=False
2004     else:
2005 caltinay 3789 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2006 gross 1841 out=False
2007     elif A.getRank() == 2:
2008     if s[0]==s[1]:
2009     for j in range(s[0]):
2010     for l in range(s[1]):
2011     if util.Lsup(A[j,l]-A[l,j])>tol:
2012 caltinay 3789 if verbose: print(("non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)))
2013 gross 1841 out=False
2014     else:
2015 caltinay 3789 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2016 gross 1841 out=False
2017     elif A.getRank() == 0:
2018     pass
2019     else:
2020 caltinay 3789 raise ValueError("Cannot check rank %s of %s."%(A.getRank(),name))
2021 gross 1841 return out
2022 jgs 148
2023 gross 1841 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
2024     """
2025 caltinay 2169 Tests two coefficients for reciprocal symmetry.
2026 jgs 148
2027 jfenwick 2625 :param name0: name of the first coefficient
2028     :type name0: ``str``
2029     :param name1: name of the second coefficient
2030     :type name1: ``str``
2031     :param verbose: if set to True or not present a report on coefficients
2032 caltinay 2169 which break the symmetry is printed
2033 jfenwick 2625 :type verbose: ``bool``
2034     :return: True if coefficients ``name0`` and ``name1`` are reciprocally
2035 caltinay 2169 symmetric.
2036 jfenwick 2625 :rtype: ``bool``
2037 gross 1841 """
2038 gross 2474 SMALL_TOLERANCE=util.EPSILON*10.
2039 gross 1841 B=self.getCoefficient(name0)
2040     C=self.getCoefficient(name1)
2041     verbose=verbose or self.__debug
2042     out=True
2043     if B.isEmpty() and not C.isEmpty():
2044 caltinay 3789 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2045 gross 1841 out=False
2046     elif not B.isEmpty() and C.isEmpty():
2047 caltinay 3789 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2048 gross 1841 out=False
2049     elif not B.isEmpty() and not C.isEmpty():
2050     sB=B.getShape()
2051     sC=C.getShape()
2052 gross 2474 tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2053 gross 1841 if len(sB) != len(sC):
2054 caltinay 3789 if verbose: print(("non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))))
2055 gross 1841 out=False
2056     else:
2057     if len(sB)==0:
2058     if util.Lsup(B-C)>tol:
2059 caltinay 3789 if verbose: print(("non-symmetric problem because %s!=%s"%(name0,name1)))
2060 gross 1841 out=False
2061     elif len(sB)==1:
2062     if sB[0]==sC[0]:
2063     for j in range(sB[0]):
2064     if util.Lsup(B[j]-C[j])>tol:
2065 caltinay 3789 if verbose: print(("non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)))
2066 gross 1841 out=False
2067     else:
2068 caltinay 3789 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2069 gross 1841 elif len(sB)==3:
2070     if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
2071     for i in range(sB[0]):
2072     for j in range(sB[1]):
2073     for k in range(sB[2]):
2074     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
2075 caltinay 3789 if verbose: print(("non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)))
2076 gross 1841 out=False
2077     else:
2078 caltinay 3789 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2079 jgs 148 else:
2080 caltinay 3789 raise ValueError("Cannot check rank %s of %s and %s."%(len(sB),name0,name1))
2081 gross 1841 return out
2082 jgs 148
2083 caltinay 2169 def getCoefficient(self,name):
2084 jgs 108 """
2085 jfenwick 2625 Returns the value of the coefficient ``name``.
2086 jgs 108
2087 jfenwick 2625 :param name: name of the coefficient requested
2088     :type name: ``string``
2089     :return: the value of the coefficient
2090     :rtype: `Data`
2091     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2092 jgs 148 """
2093     if self.hasCoefficient(name):
2094 gross 1841 return self.__COEFFICIENTS[name].getValue()
2095 jgs 148 else:
2096 caltinay 3789 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2097 jgs 108
2098 caltinay 2169 def hasCoefficient(self,name):
2099 jgs 148 """
2100 jfenwick 2625 Returns True if ``name`` is the name of a coefficient.
2101 jgs 108
2102 jfenwick 2625 :param name: name of the coefficient enquired
2103     :type name: ``string``
2104     :return: True if ``name`` is the name of a coefficient of the general PDE,
2105 caltinay 2169 False otherwise
2106 jfenwick 2625 :rtype: ``bool``
2107 jgs 148 """
2108 caltinay 3789 return name in self.__COEFFICIENTS
2109 jgs 108
2110 caltinay 2169 def createCoefficient(self, name):
2111 jgs 148 """
2112 jfenwick 2625 Creates a `Data` object corresponding to coefficient
2113     ``name``.
2114 jgs 108
2115 jfenwick 2625 :return: the coefficient ``name`` initialized to 0
2116     :rtype: `Data`
2117     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2118 jgs 148 """
2119     if self.hasCoefficient(name):
2120     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
2121     else:
2122 caltinay 3789 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2123 jgs 148
2124 caltinay 2169 def getFunctionSpaceForCoefficient(self,name):
2125 jgs 148 """
2126 jfenwick 2625 Returns the `FunctionSpace` to be used for
2127     coefficient ``name``.
2128 jgs 148
2129 jfenwick 2625 :param name: name of the coefficient enquired
2130     :type name: ``string``
2131     :return: the function space to be used for coefficient ``name``
2132     :rtype: `FunctionSpace`
2133     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2134 jgs 148 """
2135     if self.hasCoefficient(name):
2136 gross 1841 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
2137 jgs 148 else:
2138 caltinay 3789 raise ValueError("unknown coefficient %s requested"%name)
2139 gross 1841
2140 caltinay 2169 def getShapeOfCoefficient(self,name):
2141 jgs 148 """
2142 jfenwick 2625 Returns the shape of the coefficient ``name``.
2143 jgs 148
2144 jfenwick 2625 :param name: name of the coefficient enquired
2145     :type name: ``string``
2146     :return: the shape of the coefficient ``name``
2147     :rtype: ``tuple`` of ``int``
2148     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2149 jgs 148 """
2150     if self.hasCoefficient(name):
2151 gross 1841 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
2152 jgs 148 else:
2153 caltinay 3789 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2154 jgs 148
2155 caltinay 2169 def resetAllCoefficients(self):
2156 jgs 148 """
2157 caltinay 2169 Resets all coefficients to their default values.
2158 jgs 148 """
2159 caltinay 3789 for i in list(self.__COEFFICIENTS.keys()):
2160 gross 1841 self.__COEFFICIENTS[i].resetValue()
2161 jgs 148
2162 caltinay 2169 def alteredCoefficient(self,name):
2163 jgs 148 """
2164 jfenwick 2625 Announces that coefficient ``name`` has been changed.
2165 jgs 148
2166 jfenwick 2625 :param name: name of the coefficient affected
2167     :type name: ``string``
2168     :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2169     :note: if ``name`` is q or r, the method will not trigger a rebuild of the
2170 caltinay 2169 system as constraints are applied to the solved system.
2171 jgs 148 """
2172     if self.hasCoefficient(name):
2173     self.trace("Coefficient %s has been altered."%name)
2174 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
2175 gross 1841 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2176     if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2177 jgs 148 else:
2178 caltinay 3789 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2179 jgs 148
2180 caltinay 2169 def validSolution(self):
2181 gross 1841 """
2182 caltinay 2169 Marks the solution as valid.
2183 gross 1841 """
2184     self.__is_solution_valid=True
2185 jgs 149
2186 caltinay 2169 def invalidateSolution(self):
2187 gross 1841 """
2188 caltinay 2169 Indicates the PDE has to be resolved if the solution is requested.
2189 gross 1841 """
2190     self.trace("System will be resolved.")
2191     self.__is_solution_valid=False
2192 jgs 148
2193 gross 1841 def isSolutionValid(self):
2194     """
2195 caltinay 2169 Returns True if the solution is still valid.
2196 gross 1841 """
2197 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2198 gross 2474 if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2199     self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2200     self.invalidateSolution()
2201 gross 1841 return self.__is_solution_valid
2202    
2203     def validOperator(self):
2204     """
2205 caltinay 2169 Marks the operator as valid.
2206 gross 1841 """
2207     self.__is_operator_valid=True
2208    
2209 caltinay 2169 def invalidateOperator(self):
2210 gross 1841 """
2211 caltinay 2169 Indicates the operator has to be rebuilt next time it is used.
2212 gross 1841 """
2213     self.trace("Operator will be rebuilt.")
2214     self.invalidateSolution()
2215     self.__is_operator_valid=False
2216    
2217     def isOperatorValid(self):
2218     """
2219 caltinay 2169 Returns True if the operator is still valid.
2220 gross 1841 """
2221 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
2222     if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
2223 gross 1841 return self.__is_operator_valid
2224    
2225     def validRightHandSide(self):
2226     """
2227 caltinay 2169 Marks the right hand side as valid.
2228 gross 1841 """
2229     self.__is_RHS_valid=True
2230    
2231 caltinay 2169 def invalidateRightHandSide(self):
2232 gross 1841 """
2233 caltinay 2169 Indicates the right hand side has to be rebuilt next time it is used.
2234 gross 1841 """
2235 gross 2535 self.trace("Right hand side has to be rebuilt.")
2236 gross 1841 self.invalidateSolution()
2237     self.__is_RHS_valid=False
2238    
2239     def isRightHandSideValid(self):
2240     """
2241 caltinay 2169 Returns True if the operator is still valid.
2242 gross 1841 """
2243 gross 2535 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
2244 gross 1841 return self.__is_RHS_valid
2245    
2246 caltinay 2169 def invalidateSystem(self):
2247 gross 1841 """
2248 caltinay 2169 Announces that everything has to be rebuilt.
2249 gross 1841 """
2250     self.invalidateSolution()
2251     self.invalidateOperator()
2252     self.invalidateRightHandSide()
2253    
2254     def isSystemValid(self):
2255     """
2256 caltinay 2169 Returns True if the system (including solution) is still vaild.
2257 gross 1841 """
2258     return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
2259    
2260 caltinay 2169 def initializeSystem(self):
2261 gross 1841 """
2262 caltinay 2169 Resets the system clearing the operator, right hand side and solution.
2263 gross 1841 """
2264     self.trace("New System has been created.")
2265 gross 2474 self.__operator_type=None
2266 gross 2535 self.setSystemStatus()
2267 gross 1841 self.__operator=escript.Operator()
2268     self.__righthandside=escript.Data()
2269     self.__solution=escript.Data()
2270     self.invalidateSystem()
2271    
2272 caltinay 2169 def getOperator(self):
2273 gross 1841 """
2274 caltinay 2169 Returns the operator of the linear problem.
2275 gross 1841
2276 jfenwick 2625 :return: the operator of the problem
2277 gross 1841 """
2278     return self.getSystem()[0]
2279    
2280 caltinay 2169 def getRightHandSide(self):
2281 gross 1841 """
2282 caltinay 2169 Returns the right hand side of the linear problem.
2283 gross 1841
2284 jfenwick 2625 :return: the right hand side of the problem
2285     :rtype: `Data`
2286 gross 1841 """
2287     return self.getSystem()[1]
2288    
2289 caltinay 2169 def createRightHandSide(self):
2290 gross 1841 """
2291 caltinay 2169 Returns an instance of a new right hand side.
2292 gross 1841 """
2293     self.trace("New right hand side is allocated.")
2294     if self.getNumEquations()>1:
2295     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2296     else:
2297     return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2298    
2299 caltinay 2169 def createSolution(self):
2300 gross 1841 """
2301 caltinay 2169 Returns an instance of a new solution.
2302 gross 1841 """
2303     self.trace("New solution is allocated.")
2304     if self.getNumSolutions()>1:
2305     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2306     else:
2307     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
2308    
2309 caltinay 2169 def resetSolution(self):
2310 gross 1841 """
2311 caltinay 2169 Sets the solution to zero.
2312 gross 1841 """
2313     if self.__solution.isEmpty():
2314     self.__solution=self.createSolution()
2315     else:
2316     self.__solution.setToZero()
2317     self.trace("Solution is reset to zero.")
2318 caltinay 2169
2319 gross 2987 def setSolution(self,u, validate=True):
2320 gross 1841 """
2321 gross 2474 Sets the solution assuming that makes the system valid with the tolrance
2322     defined by the solver options
2323 gross 1841 """
2324 gross 2987 if validate:
2325     self.__solution_rtol=self.getSolverOptions().getTolerance()
2326     self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2327     self.validSolution()
2328 gross 1841 self.__solution=u
2329     def getCurrentSolution(self):
2330     """
2331 caltinay 2169 Returns the solution in its current state.
2332 gross 1841 """
2333     if self.__solution.isEmpty(): self.__solution=self.createSolution()
2334     return self.__solution
2335    
2336 caltinay 2169 def resetRightHandSide(self):
2337 gross 1841 """
2338 caltinay 2169 Sets the right hand side to zero.
2339 gross 1841 """
2340     if self.__righthandside.isEmpty():
2341     self.__righthandside=self.createRightHandSide()
2342     else:
2343     self.__righthandside.setToZero()
2344     self.trace("Right hand side is reset to zero.")
2345    
2346     def getCurrentRightHandSide(self):
2347     """
2348 caltinay 2169 Returns the right hand side in its current state.
2349 gross 1841 """
2350     if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2351     return self.__righthandside
2352    
2353 caltinay 2169 def resetOperator(self):
2354 gross 1841 """
2355 caltinay 2169 Makes sure that the operator is instantiated and returns it initialized
2356     with zeros.
2357 gross 1841 """
2358 gross 2474 if self.getOperatorType() == None:
2359 gross 1841 if self.isUsingLumping():
2360     self.__operator=self.createSolution()
2361     else:
2362     self.__operator=self.createOperator()
2363 gross 2474 self.__operator_type=self.getRequiredOperatorType()
2364 gross 1841 else:
2365     if self.isUsingLumping():
2366     self.__operator.setToZero()
2367     else: