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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3815 - (show annotations)
Thu Feb 9 00:27:46 2012 UTC (7 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 168337 byte(s)
Merging trunk 3814 into symbolic to get ripley.

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