/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


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