/[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 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 167699 byte(s)
Fast forward to latest trunk revision 3788.

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