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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4622 - (show annotations)
Fri Jan 17 04:55:41 2014 UTC (5 years, 3 months ago) by sshaw
File MIME type: text/x-python
File size: 176475 byte(s)
Added dirac support to ripley, added interface for custom assemblers for ripleydomains (also added the custom assembler for 2D VTI waves), changed synthetic_VTI example to use the new, faster custom assembler

1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2013 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-2013 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 from . import escriptcpp as escore
42 from . import util
43 import math
44 import numpy
45
46 __author__="Lutz Gross, l.gross@uq.edu.au"
47
48
49 class SolverOptions(object):
50 """
51 this class defines the solver options for a linear or non-linear solver.
52
53 The option also supports the handling of diagnostic informations.
54
55 Typical usage is
56
57 ::
58
59 opts=SolverOptions()
60 print(opts)
61 opts.resetDiagnostics()
62 u=solver(opts)
63 print("number of iteration steps: =",opts.getDiagnostics("num_iter"))
64
65 :cvar DEFAULT: The default method used to solve the system of linear equations
66 :cvar DIRECT: The direct solver based on LDU factorization
67 :cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
68 :cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
69 :cvar CR: The conjugate residual method
70 :cvar CGS: The conjugate gradient square method
71 :cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
72 :cvar TFQMR: Transpose Free Quasi Minimal Residual method
73 :cvar MINRES: Minimum residual method
74 :cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
75 :cvar ILUT: The incomplete LU factorization preconditioner with fill-in
76 :cvar JACOBI: The Jacobi preconditioner
77 :cvar GMRES: The Gram-Schmidt minimum residual method
78 :cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
79 :cvar ROWSUM_LUMPING: Matrix lumping using row sum
80 :cvar HRZ_LUMPING: Matrix lumping using the HRZ approach
81 :cvar NO_REORDERING: No matrix reordering allowed
82 :cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
83 :cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
84 :cvar PASO: PASO solver package
85 :cvar SCSL: SGI SCSL solver library
86 :cvar MKL: Intel's MKL solver library
87 :cvar UMFPACK: The UMFPACK library
88 :cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
89 :cvar ITERATIVE: The default iterative solver
90 :cvar AMG: Algebraic Multi Grid
91 :cvar AMLI: Algebraic Multi Level Iteration
92 :cvar REC_ILU: recursive ILU0
93 :cvar RILU: relaxed ILU0
94 :cvar GAUSS_SEIDEL: Gauss-Seidel preconditioner
95 :cvar DEFAULT_REORDERING: the reordering method recommended by the solver
96 :cvar SUPER_LU: the Super_LU solver package
97 :cvar PASTIX: the Pastix direct solver_package
98 :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
99 :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
100 :cvar AGGREGATION_COARSENING: AMG and AMLI coarsening using (symmetric) aggregation
101 :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
102 :cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
103 :cvar NO_PRECONDITIONER: no preconditioner is applied.
104 :cvar CLASSIC_INTERPOLATION_WITH_FF_COUPLING: classical interpolation in AMG with enforced
105 :cvar CLASSIC_INTERPOLATION: classical interpolation in AMG
106 :cvar DIRECT_INTERPOLATION: direct interploation in AMG
107 :cvar BOOMERAMG: Boomer AMG in hypre library
108 :cvar CIJP_FIXED_RANDOM_COARSENING: BoomerAMG parallel coarsening method CIJP by using fixed random vector
109 :cvar CIJP_COARSENING: BoomerAMG parallel coarsening method CIJP
110 :cvar PASO_FALGOUT_COARSENING: BoomerAMG parallel coarsening method falgout
111 :cvar PASO_PMIS_COARSENING: BoomerAMG parallel coarsening method PMIS
112 :cvar PASO_HMIS_COARSENING: BoomerAMG parallel coarsening method HMIS
113 :cvar BACKWARD_EULER: backward Euler scheme
114 :cvar CRANK_NICOLSON: Crank-Nicolson scheme
115 :cvar LINEAR_CRANK_NICOLSON: linerized Crank-Nicolson scheme
116
117 """
118 DEFAULT= 0
119 DIRECT= 1
120 CHOLEVSKY= 2
121 PCG= 3
122 CR= 4
123 CGS= 5
124 BICGSTAB= 6
125 ILU0= 8
126 ILUT= 9
127 JACOBI= 10
128 GMRES= 11
129 PRES20= 12
130 LUMPING= 13
131 ROWSUM_LUMPING= 13
132 HRZ_LUMPING= 14
133 NO_REORDERING= 17
134 MINIMUM_FILL_IN= 18
135 NESTED_DISSECTION= 19
136 MKL= 15
137 UMFPACK= 16
138 ITERATIVE= 20
139 PASO= 21
140 AMG= 22
141 REC_ILU = 23
142 TRILINOS = 24
143 NONLINEAR_GMRES = 25
144 TFQMR = 26
145 MINRES = 27
146 GAUSS_SEIDEL=28
147 RILU=29
148 DEFAULT_REORDERING=30
149 SUPER_LU=31
150 PASTIX=32
151 YAIR_SHAPIRA_COARSENING=33
152 RUGE_STUEBEN_COARSENING=34
153 AGGREGATION_COARSENING=35
154 NO_PRECONDITIONER=36
155 MIN_COARSE_MATRIX_SIZE=37
156 AMLI=38
157 STANDARD_COARSENING=39
158 CLASSIC_INTERPOLATION_WITH_FF_COUPLING=50
159 CLASSIC_INTERPOLATION=51
160 DIRECT_INTERPOLATION=52
161 BOOMERAMG=60
162 CIJP_FIXED_RANDOM_COARSENING=61
163 CIJP_COARSENING=62
164 FALGOUT_COARSENING=63
165 PMIS_COARSENING=64
166 HMIS_COARSENING=65
167 LINEAR_CRANK_NICOLSON=66
168 CRANK_NICOLSON=67
169 BACKWARD_EULER=68
170
171 def __init__(self):
172 self.setLevelMax()
173 self.setCoarseningThreshold()
174 self.setSmoother()
175 self.setNumSweeps()
176 self.setNumPreSweeps()
177 self.setNumPostSweeps()
178 self.setTolerance()
179 self.setAbsoluteTolerance()
180 self.setInnerTolerance()
181 self.setDropTolerance()
182 self.setDropStorage()
183 self.setIterMax()
184 self.setInnerIterMax()
185 self.setTruncation()
186 self.setRestart()
187 self.setSymmetry()
188 self.setVerbosity()
189 self.setInnerToleranceAdaption()
190 self.setAcceptanceConvergenceFailure()
191 self.setReordering()
192 self.setPackage()
193 self.setSolverMethod()
194 self.setPreconditioner()
195 self.setCoarsening()
196 self.setMinCoarseMatrixSize()
197 self.setRelaxationFactor()
198 self.setLocalPreconditionerOff()
199 self.resetDiagnostics(all=True)
200 self.setMinCoarseMatrixSparsity()
201 self.setNumRefinements()
202 self.setNumCoarseMatrixRefinements()
203 self.setUsePanel()
204 self.setDiagonalDominanceThreshold()
205 self.setAMGInterpolation()
206 self.setCycleType()
207 self.setODESolver()
208
209
210 def __str__(self):
211 return self.getSummary()
212 def getSummary(self):
213 """
214 Returns a string reporting the current settings
215 """
216 out="Solver Package: %s"%(self.getName(self.getPackage()))
217 out+="\nVerbosity = %s"%self.isVerbose()
218 out+="\nAccept failed convergence = %s"%self.acceptConvergenceFailure()
219 out+="\nRelative tolerance = %e"%self.getTolerance()
220 out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
221 out+="\nSymmetric problem = %s"%self.isSymmetric()
222 out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
223 # out+="\nInner tolerance = %e"%self.getInnerTolerance()
224 # out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
225
226 if self.getPackage() in (self.DEFAULT, self.PASO):
227 out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
228 if self.getSolverMethod() == self.GMRES:
229 out+="\nTruncation = %s"%self.getTruncation()
230 out+="\nRestart = %s"%self.getRestart()
231 if self.getSolverMethod() == self.AMG:
232 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
233 out+="\nMaximum number of levels = %s"%self.getLevelMax()
234 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
235 out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
236 out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
237 out+="\nApply preconditioner locally = %s"%self.useLocalPreconditioner()
238 if self.getPreconditioner() == self.AMG:
239 out+="\nMaximum number of levels = %s"%self.getLevelMax()
240 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
241 out+="\nMinimal sparsity on coarsest level = %e"%self.getMinCoarseMatrixSparsity()
242 out+="\nSmoother = %s"%self.getName(self.getSmoother())
243 out+="\nMinimum size of the coarsest level matrix = %e"%self.getMinCoarseMatrixSize()
244 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
245 out+="\nNumber of refinement steps in coarsest level solver = %d"%self.getNumCoarseMatrixRefinements()
246 out+="\nUse node panel = %s"%self.usePanel()
247 out+="\nInterpolation = %s"%(self.getName(self.getAMGInterpolation()))
248 out+="\nThreshold for diagonal dominant rows = %s"%self.getDiagonalDominanceThreshold()
249
250 if self.getPreconditioner() == self.AMLI:
251 out+="\nMaximum number of levels = %s"%self.getLevelMax()
252 out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
253 out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
254 out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
255 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
256 if self.getPreconditioner() == self.BOOMERAMG:
257 out+="\nMaximum number of levels = %s"%self.getLevelMax()
258 out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
259 out+="\nThreshold for diagonal dominant rows = %s"%setDiagonalDominanceThreshold()
260 out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
261 out+="\nV-cycle (1) or W-cyle (2) = %s"%self.getCycleType()
262 out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
263 out+="\nSmoother = %s"%self.getName(self.getSmoother())
264 if self.getPreconditioner() == self.GAUSS_SEIDEL:
265 out+="\nNumber of sweeps = %s"%self.getNumSweeps()
266 if self.getPreconditioner() == self.ILUT:
267 out+="\nDrop tolerance = %e"%self.getDropTolerance()
268 out+="\nStorage increase = %e"%self.getDropStorage()
269 if self.getPreconditioner() == self.RILU:
270 out+="\nRelaxation factor = %e"%self.getRelaxationFactor()
271 out+="\nODE solver = %s"%self.getName(self.getODESolver())
272 return out
273
274 def getName(self,key):
275 """
276 returns the name of a given key
277
278 :param key: a valid key
279 """
280 if key == self.DEFAULT: return "DEFAULT"
281 if key == self.DIRECT: return "DIRECT"
282 if key == self.CHOLEVSKY: return "CHOLEVSKY"
283 if key == self.PCG: return "PCG"
284 if key == self.CR: return "CR"
285 if key == self.CGS: return "CGS"
286 if key == self.BICGSTAB: return "BICGSTAB"
287 if key == self.ILU0: return "ILU0:"
288 if key == self.ILUT: return "ILUT"
289 if key == self.JACOBI: return "JACOBI"
290 if key == self.GMRES: return "GMRES"
291 if key == self.PRES20: return "PRES20"
292 if key == self.ROWSUM_LUMPING: return "ROWSUM_LUMPING"
293 if key == self.HRZ_LUMPING: return "HRZ_LUMPING"
294 if key == self.NO_REORDERING: return "NO_REORDERING"
295 if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
296 if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
297 if key == self.MKL: return "MKL"
298 if key == self.UMFPACK: return "UMFPACK"
299 if key == self.ITERATIVE: return "ITERATIVE"
300 if key == self.PASO: return "PASO"
301 if key == self.AMG: return "AMG"
302 if key == self.AMLI: return "AMLI"
303 if key == self.REC_ILU: return "REC_ILU"
304 if key == self.TRILINOS: return "TRILINOS"
305 if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
306 if key == self.TFQMR: return "TFQMR"
307 if key == self.MINRES: return "MINRES"
308 if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
309 if key == self.RILU: return "RILU"
310 if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
311 if key == self.SUPER_LU: return "SUPER_LU"
312 if key == self.PASTIX: return "PASTIX"
313 if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
314 if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
315 if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
316 if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
317 if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
318 if key == self.CLASSIC_INTERPOLATION_WITH_FF_COUPLING: return "CLASSIC_INTERPOLATION_WITH_FF"
319 if key == self.CLASSIC_INTERPOLATION: return "CLASSIC_INTERPOLATION"
320 if key == self.DIRECT_INTERPOLATION: return "DIRECT_INTERPOLATION"
321 if key == self.BOOMERAMG: return "BOOMERAMG"
322 if key == self.CIJP_FIXED_RANDOM_COARSENING: return "CIJP_FIXED_RANDOM_COARSENING"
323 if key == self.CIJP_COARSENING: return "CIJP_COARSENING"
324 if key == self.FALGOUT_COARSENING: return "FALGOUT_COARSENING"
325 if key == self.PMIS_COARSENING: return "PMIS_COARSENING"
326 if key == self.HMIS_COARSENING: return "HMIS_COARSENING"
327 if key == self.LINEAR_CRANK_NICOLSON: return "LINEAR_CRANK_NICOLSON"
328 if key == self.CRANK_NICOLSON: return "CRANK_NICOLSON"
329 if key == self.BACKWARD_EULER: return "BACKWARD_EULER"
330
331
332
333 def resetDiagnostics(self,all=False):
334 """
335 resets the diagnostics
336
337 :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
338 :type all: ``bool``
339 """
340 self.__num_iter=None
341 self.__num_level=None
342 self.__num_inner_iter=None
343 self.__time=None
344 self.__set_up_time=None
345 self.__net_time=None
346 self.__residual_norm=None
347 self.__converged=None
348 self.__preconditioner_size=-1
349 self.__time_step_backtracking_used=None
350 if all:
351 self.__cum_num_inner_iter=0
352 self.__cum_num_iter=0
353 self.__cum_time=0
354 self.__cum_set_up_time=0
355 self.__cum_net_time=0
356
357 def _updateDiagnostics(self, name, value):
358 """
359 Updates diagnostic information
360
361 :param name: name of diagnostic information
362 :type name: ``str`` in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
363 :param value: new value of the diagnostic information
364 :note: this function is used by a solver to report diagnostics informations.
365 """
366 if name == "num_iter":
367 self.__num_iter=int(value)
368 self.__cum_num_iter+=self.__num_iter
369 if name == "num_level":
370 self.__num_level=int(value)
371 if name == "num_inner_iter":
372 self.__num_inner_iter=int(value)
373 self.__cum_num_inner_iter+=self.__num_inner_iter
374 if name == "time":
375 self.__time=float(value)
376 self.__cum_time+=self.__time
377 if name == "set_up_time":
378 self.__set_up_time=float(value)
379 self.__cum_set_up_time+=self.__set_up_time
380 if name == "net_time":
381 self.__net_time=float(value)
382 self.__cum_net_time+=self.__net_time
383 if name == "residual_norm":
384 self.__residual_norm=float(value)
385 if name == "converged":
386 self.__converged = (value == True)
387 if name == "time_step_backtracking_used":
388 self.__time_step_backtracking_used = (value == True)
389 if name == "coarse_level_sparsity":
390 self.__coarse_level_sparsity = value
391 if name == "num_coarse_unknowns":
392 self.__num_coarse_unknowns= value
393
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("minimum size of the coarsest level matrix must be non-negative.")
477 self.__MinCoarseMatrixSize=size
478
479 def getMinCoarseMatrixSize(self):
480 """
481 Returns the minimum 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 coefficient alters the operator of
1305 the PDE
1306 :cvar RIGHTHANDSIDE: indicator that the coefficient alters the right
1307 hand side of the PDE
1308 :cvar BOTH: indicator that 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 def resetRightHandSideCoefficients(self):
1707 """
1708 Resets all coefficients defining the right hand side
1709 """
1710 for name in self.__COEFFICIENTS:
1711 if self.__COEFFICIENTS[name].altering == PDECoef.RIGHTHANDSIDE :
1712 self.__COEFFICIENTS[name].resetValue()
1713 self.trace("coefficient %s has been reset."%name)
1714
1715 def getDomain(self):
1716 """
1717 Returns the domain of the PDE.
1718
1719 :return: the domain of the PDE
1720 :rtype: `Domain`
1721 """
1722 return self.__domain
1723 def getDomainStatus(self):
1724 """
1725 Return the status indicator of the domain
1726 """
1727 return self.getDomain().getStatus()
1728
1729 def getSystemStatus(self):
1730 """
1731 Return the domain status used to build the current system
1732 """
1733 return self.__system_status
1734 def setSystemStatus(self,status=None):
1735 """
1736 Sets the system status to ``status`` if ``status`` is not present the
1737 current status of the domain is used.
1738 """
1739 if status == None:
1740 self.__system_status=self.getDomainStatus()
1741 else:
1742 self.__system_status=status
1743
1744 def getDim(self):
1745 """
1746 Returns the spatial dimension of the PDE.
1747
1748 :return: the spatial dimension of the PDE domain
1749 :rtype: ``int``
1750 """
1751 return self.getDomain().getDim()
1752
1753 def getNumEquations(self):
1754 """
1755 Returns the number of equations.
1756
1757 :return: the number of equations
1758 :rtype: ``int``
1759 :raise UndefinedPDEError: if the number of equations is not specified yet
1760 """
1761 if self.__numEquations==None:
1762 if self.__numSolutions==None:
1763 raise UndefinedPDEError("Number of equations is undefined. Please specify argument numEquations.")
1764 else:
1765 self.__numEquations=self.__numSolutions
1766 return self.__numEquations
1767
1768 def getNumSolutions(self):
1769 """
1770 Returns the number of unknowns.
1771
1772 :return: the number of unknowns
1773 :rtype: ``int``
1774 :raise UndefinedPDEError: if the number of unknowns is not specified yet
1775 """
1776 if self.__numSolutions==None:
1777 if self.__numEquations==None:
1778 raise UndefinedPDEError("Number of solution is undefined. Please specify argument numSolutions.")
1779 else:
1780 self.__numSolutions=self.__numEquations
1781 return self.__numSolutions
1782
1783 def reduceEquationOrder(self):
1784 """
1785 Returns the status of order reduction for the equation.
1786
1787 :return: True if reduced interpolation order is used for the
1788 representation of the equation, False otherwise
1789 :rtype: `bool`
1790 """
1791 return self.__reduce_equation_order
1792
1793 def reduceSolutionOrder(self):
1794 """
1795 Returns the status of order reduction for the solution.
1796
1797 :return: True if reduced interpolation order is used for the
1798 representation of the solution, False otherwise
1799 :rtype: `bool`
1800 """
1801 return self.__reduce_solution_order
1802
1803 def getFunctionSpaceForEquation(self):
1804 """
1805 Returns the `FunctionSpace` used to discretize
1806 the equation.
1807
1808 :return: representation space of equation
1809 :rtype: `FunctionSpace`
1810 """
1811 if self.reduceEquationOrder():
1812 return escore.ReducedSolution(self.getDomain())
1813 else:
1814 return escore.Solution(self.getDomain())
1815
1816 def getFunctionSpaceForSolution(self):
1817 """
1818 Returns the `FunctionSpace` used to represent
1819 the solution.
1820
1821 :return: representation space of solution
1822 :rtype: `FunctionSpace`
1823 """
1824 if self.reduceSolutionOrder():
1825 return escore.ReducedSolution(self.getDomain())
1826 else:
1827 return escore.Solution(self.getDomain())
1828
1829 # ==========================================================================
1830 # solver settings:
1831 # ==========================================================================
1832 def setSolverOptions(self,options=None):
1833 """
1834 Sets the solver options.
1835
1836 :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1837 :type options: `SolverOptions` or ``None``
1838 :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1839 """
1840 if options==None:
1841 self.__solver_options=SolverOptions()
1842 elif isinstance(options, SolverOptions):
1843 self.__solver_options=options
1844 else:
1845 raise ValueError("options must be a SolverOptions object.")
1846 self.__solver_options.setSymmetry(self.__sym)
1847
1848 def getSolverOptions(self):
1849 """
1850 Returns the solver options
1851
1852 :rtype: `SolverOptions`
1853 """
1854 self.__solver_options.setSymmetry(self.__sym)
1855 return self.__solver_options
1856
1857 def isUsingLumping(self):
1858 """
1859 Checks if matrix lumping is the current solver method.
1860
1861 :return: True if the current solver method is lumping
1862 :rtype: ``bool``
1863 """
1864 return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
1865 # ==========================================================================
1866 # symmetry flag:
1867 # ==========================================================================
1868 def isSymmetric(self):
1869 """
1870 Checks if symmetry is indicated.
1871
1872 :return: True if a symmetric PDE is indicated, False otherwise
1873 :rtype: ``bool``
1874 :note: the method is equivalent to use getSolverOptions().isSymmetric()
1875 """
1876 self.getSolverOptions().isSymmetric()
1877
1878 def setSymmetryOn(self):
1879 """
1880 Sets the symmetry flag.
1881 :note: The method overwrites the symmetry flag set by the solver options
1882 """
1883 self.__sym=True
1884 self.getSolverOptions().setSymmetryOn()
1885
1886 def setSymmetryOff(self):
1887 """
1888 Clears the symmetry flag.
1889 :note: The method overwrites the symmetry flag set by the solver options
1890 """
1891 self.__sym=False
1892 self.getSolverOptions().setSymmetryOff()
1893
1894 def setSymmetry(self,flag=False):
1895 """
1896 Sets the symmetry flag to ``flag``.
1897
1898 :param flag: If True, the symmetry flag is set otherwise reset.
1899 :type flag: ``bool``
1900 :note: The method overwrites the symmetry flag set by the solver options
1901 """
1902 self.getSolverOptions().setSymmetry(flag)
1903 # ==========================================================================
1904 # function space handling for the equation as well as the solution
1905 # ==========================================================================
1906 def setReducedOrderOn(self):
1907 """
1908 Switches reduced order on for solution and equation representation.
1909
1910 :raise RuntimeError: if order reduction is altered after a coefficient has
1911 been set
1912 """
1913 self.setReducedOrderForSolutionOn()
1914 self.setReducedOrderForEquationOn()
1915
1916 def setReducedOrderOff(self):
1917 """
1918 Switches reduced order off for solution and equation representation
1919
1920 :raise RuntimeError: if order reduction is altered after a coefficient has
1921 been set
1922 """
1923 self.setReducedOrderForSolutionOff()
1924 self.setReducedOrderForEquationOff()
1925
1926 def setReducedOrderTo(self,flag=False):
1927 """
1928 Sets order reduction state for both solution and equation representation
1929 according to flag.
1930
1931 :param flag: if True, the order reduction is switched on for both solution
1932 and equation representation, otherwise or if flag is not
1933 present order reduction is switched off
1934 :type flag: ``bool``
1935 :raise RuntimeError: if order reduction is altered after a coefficient has
1936 been set
1937 """
1938 self.setReducedOrderForSolutionTo(flag)
1939 self.setReducedOrderForEquationTo(flag)
1940
1941
1942 def setReducedOrderForSolutionOn(self):
1943 """
1944 Switches reduced order on for solution representation.
1945
1946 :raise RuntimeError: if order reduction is altered after a coefficient has
1947 been set
1948 """
1949 if not self.__reduce_solution_order:
1950 if self.__altered_coefficients:
1951 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1952 self.trace("Reduced order is used for solution representation.")
1953 self.__reduce_solution_order=True
1954 self.initializeSystem()
1955
1956 def setReducedOrderForSolutionOff(self):
1957 """
1958 Switches reduced order off for solution representation
1959
1960 :raise RuntimeError: if order reduction is altered after a coefficient has
1961 been set.
1962 """
1963 if self.__reduce_solution_order:
1964 if self.__altered_coefficients:
1965 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1966 self.trace("Full order is used to interpolate solution.")
1967 self.__reduce_solution_order=False
1968 self.initializeSystem()
1969
1970 def setReducedOrderForSolutionTo(self,flag=False):
1971 """
1972 Sets order reduction state for solution representation according to flag.
1973
1974 :param flag: if flag is True, the order reduction is switched on for
1975 solution representation, otherwise or if flag is not present
1976 order reduction is switched off
1977 :type flag: ``bool``
1978 :raise RuntimeError: if order reduction is altered after a coefficient has
1979 been set
1980 """
1981 if flag:
1982 self.setReducedOrderForSolutionOn()
1983 else:
1984 self.setReducedOrderForSolutionOff()
1985
1986 def setReducedOrderForEquationOn(self):
1987 """
1988 Switches reduced order on for equation representation.
1989
1990 :raise RuntimeError: if order reduction is altered after a coefficient has
1991 been set
1992 """
1993 if not self.__reduce_equation_order:
1994 if self.__altered_coefficients:
1995 raise RuntimeError("order cannot be altered after coefficients have been defined.")
1996 self.trace("Reduced order is used for test functions.")
1997 self.__reduce_equation_order=True
1998 self.initializeSystem()
1999
2000 def setReducedOrderForEquationOff(self):
2001 """
2002 Switches reduced order off for equation representation.
2003
2004 :raise RuntimeError: if order reduction is altered after a coefficient has
2005 been set
2006 """
2007 if self.__reduce_equation_order:
2008 if self.__altered_coefficients:
2009 raise RuntimeError("order cannot be altered after coefficients have been defined.")
2010 self.trace("Full order is used for test functions.")
2011 self.__reduce_equation_order=False
2012 self.initializeSystem()
2013
2014 def setReducedOrderForEquationTo(self,flag=False):
2015 """
2016 Sets order reduction state for equation representation according to flag.
2017
2018 :param flag: if flag is True, the order reduction is switched on for
2019 equation representation, otherwise or if flag is not present
2020 order reduction is switched off
2021 :type flag: ``bool``
2022 :raise RuntimeError: if order reduction is altered after a coefficient has
2023 been set
2024 """
2025 if flag:
2026 self.setReducedOrderForEquationOn()
2027 else:
2028 self.setReducedOrderForEquationOff()
2029 def getOperatorType(self):
2030 """
2031 Returns the current system type.
2032 """
2033 return self.__operator_type
2034
2035 def checkSymmetricTensor(self,name,verbose=True):
2036 """
2037 Tests a coefficient for symmetry.
2038
2039 :param name: name of the coefficient
2040 :type name: ``str``
2041 :param verbose: if set to True or not present a report on coefficients
2042 which break the symmetry is printed.
2043 :type verbose: ``bool``
2044 :return: True if coefficient ``name`` is symmetric
2045 :rtype: ``bool``
2046 """
2047 SMALL_TOLERANCE=util.EPSILON*10.
2048 A=self.getCoefficient(name)
2049 verbose=verbose or self.__debug
2050 out=True
2051 if not A.isEmpty():
2052 tol=util.Lsup(A)*SMALL_TOLERANCE
2053 s=A.getShape()
2054 if A.getRank() == 4:
2055 if s[0]==s[2] and s[1] == s[3]:
2056 for i in range(s[0]):
2057 for j in range(s[1]):
2058 for k in range(s[2]):
2059 for l in range(s[3]):
2060 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
2061 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)))
2062 out=False
2063 else:
2064 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2065 out=False
2066 elif A.getRank() == 2:
2067 if s[0]==s[1]:
2068 for j in range(s[0]):
2069 for l in range(s[1]):
2070 if util.Lsup(A[j,l]-A[l,j])>tol:
2071 if verbose: print(("non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)))
2072 out=False
2073 else:
2074 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
2075 out=False
2076 elif A.getRank() == 0:
2077 pass
2078 else:
2079 raise ValueError("Cannot check rank %s of %s."%(A.getRank(),name))
2080 return out
2081
2082 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
2083 """
2084 Tests two coefficients for reciprocal symmetry.
2085
2086 :param name0: name of the first coefficient
2087 :type name0: ``str``
2088 :param name1: name of the second coefficient
2089 :type name1: ``str``
2090 :param verbose: if set to True or not present a report on coefficients
2091 which break the symmetry is printed
2092 :type verbose: ``bool``
2093 :return: True if coefficients ``name0`` and ``name1`` are reciprocally
2094 symmetric.
2095 :rtype: ``bool``
2096 """
2097 SMALL_TOLERANCE=util.EPSILON*10.
2098 B=self.getCoefficient(name0)
2099 C=self.getCoefficient(name1)
2100 verbose=verbose or self.__debug
2101 out=True
2102 if B.isEmpty() and not C.isEmpty():
2103 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2104 out=False
2105 elif not B.isEmpty() and C.isEmpty():
2106 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
2107 out=False
2108 elif not B.isEmpty() and not C.isEmpty():
2109 sB=B.getShape()
2110 sC=C.getShape()
2111 tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2112 if len(sB) != len(sC):
2113 if verbose: print(("non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))))
2114 out=False
2115 else:
2116 if len(sB)==0:
2117 if util.Lsup(B-C)>tol:
2118 if verbose: print(("non-symmetric problem because %s!=%s"%(name0,name1)))
2119 out=False
2120 elif len(sB)==1:
2121 if sB[0]==sC[0]:
2122 for j in range(sB[0]):
2123 if util.Lsup(B[j]-C[j])>tol:
2124 if verbose: print(("non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)))
2125 out=False
2126 else:
2127 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2128 elif len(sB)==3:
2129 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
2130 for i in range(sB[0]):
2131 for j in range(sB[1]):
2132 for k in range(sB[2]):
2133 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
2134 if verbose: print(("non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)))
2135 out=False
2136 else:
2137 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
2138 else:
2139 raise ValueError("Cannot check rank %s of %s and %s."%(len(sB),name0,name1))
2140 return out
2141
2142 def getCoefficient(self,name):
2143 """
2144 Returns the value of the coefficient ``name``.
2145
2146 :param name: name of the coefficient requested
2147 :type name: ``string``
2148 :return: the value of the coefficient
2149 :rtype: `Data`
2150 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2151 """
2152 if self.hasCoefficient(name):
2153 return self.__COEFFICIENTS[name].getValue()
2154 else:
2155 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2156
2157 def hasCoefficient(self,name):
2158 """
2159 Returns True if ``name`` is the name of a coefficient.
2160
2161 :param name: name of the coefficient enquired
2162 :type name: ``string``
2163 :return: True if ``name`` is the name of a coefficient of the general PDE,
2164 False otherwise
2165 :rtype: ``bool``
2166 """
2167 return name in self.__COEFFICIENTS
2168
2169 def createCoefficient(self, name):
2170 """
2171 Creates a `Data` object corresponding to coefficient
2172 ``name``.
2173
2174 :return: the coefficient ``name`` initialized to 0
2175 :rtype: `Data`
2176 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2177 """
2178 if self.hasCoefficient(name):
2179 return escore.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
2180 else:
2181 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2182
2183 def getFunctionSpaceForCoefficient(self,name):
2184 """
2185 Returns the `FunctionSpace` to be used for
2186 coefficient ``name``.
2187
2188 :param name: name of the coefficient enquired
2189 :type name: ``string``
2190 :return: the function space to be used for coefficient ``name``
2191 :rtype: `FunctionSpace`
2192 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2193 """
2194 if self.hasCoefficient(name):
2195 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
2196 else:
2197 raise ValueError("unknown coefficient %s requested"%name)
2198
2199 def getShapeOfCoefficient(self,name):
2200 """
2201 Returns the shape of the coefficient ``name``.
2202
2203 :param name: name of the coefficient enquired
2204 :type name: ``string``
2205 :return: the shape of the coefficient ``name``
2206 :rtype: ``tuple`` of ``int``
2207 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2208 """
2209 if self.hasCoefficient(name):
2210 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
2211 else:
2212 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2213
2214 def resetAllCoefficients(self):
2215 """
2216 Resets all coefficients to their default values.
2217 """
2218 for i in list(self.__COEFFICIENTS.keys()):
2219 self.__COEFFICIENTS[i].resetValue()
2220
2221 def alteredCoefficient(self,name):
2222 """
2223 Announces that coefficient ``name`` has been changed.
2224
2225 :param name: name of the coefficient affected
2226 :type name: ``string``
2227 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2228 :note: if ``name`` is q or r, the method will not trigger a rebuild of the
2229 system as constraints are applied to the solved system.
2230 """
2231 if self.hasCoefficient(name):
2232 self.trace("Coefficient %s has been altered."%name)
2233 if not ((name=="q" or name=="r") and self.isUsingLumping()):
2234 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2235 if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2236 else:
2237 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
2238
2239 def validSolution(self):
2240 """
2241 Marks the solution as valid.
2242 """
2243 self.__is_solution_valid=True
2244
2245 def invalidateSolution(self):
2246 """
2247 Indicates the PDE has to be resolved if the solution is requested.
2248 """
2249 self.trace("System will be resolved.")
2250 self.__is_solution_valid=False
2251
2252 def isSolutionValid(self):
2253 """
2254 Returns True if the solution is still valid.
2255 """
2256 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2257 if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2258 self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2259 self.invalidateSolution()
2260 return self.__is_solution_valid
2261
2262 def validOperator(self):
2263 """
2264 Marks the operator as valid.
2265 """
2266 self.__is_operator_valid=True
2267
2268 def invalidateOperator(self):
2269 """
2270 Indicates the operator has to be rebuilt next time it is used.
2271 """
2272 self.trace("Operator will be rebuilt.")
2273 self.invalidateSolution()
2274 self.__is_operator_valid=False
2275
2276 def isOperatorValid(self):
2277 """
2278 Returns True if the operator is still valid.
2279 """
2280 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
2281 if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
2282 return self.__is_operator_valid
2283
2284 def validRightHandSide(self):
2285 """
2286 Marks the right hand side as valid.
2287 """
2288 self.__is_RHS_valid=True
2289
2290 def invalidateRightHandSide(self):
2291 """
2292 Indicates the right hand side has to be rebuilt next time it is used.
2293 """
2294 self.trace("Right hand side has to be rebuilt.")
2295 self.invalidateSolution()
2296 self.__is_RHS_valid=False
2297
2298 def isRightHandSideValid(self):
2299 """
2300 Returns True if the operator is still valid.
2301 """
2302 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
2303 return self.__is_RHS_valid
2304
2305 def invalidateSystem(self):
2306 """
2307 Announces that everything has to be rebuilt.
2308 """
2309 self.invalidateSolution()
2310 self.invalidateOperator()
2311 self.invalidateRightHandSide()
2312
2313 def isSystemValid(self):
2314 """
2315 Returns True if the system (including solution) is still vaild.
2316 """
2317 return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
2318
2319 def initializeSystem(self):
2320 """
2321 Resets the system clearing the operator, right hand side and solution.
2322 """
2323 self.trace("New System has been created.")
2324 self.__operator_type=None
2325 self.setSystemStatus()
2326 self.__operator=escore.Operator()
2327 self.__righthandside=escore.Data()
2328 self.__solution=escore.Data()
2329 self.invalidateSystem()
2330
2331 def getOperator(self):
2332 """
2333 Returns the operator of the linear problem.
2334
2335 :return: the operator of the problem
2336 """
2337 return self.getSystem()[0]
2338
2339 def getRightHandSide(self):
2340 """
2341 Returns the right hand side of the linear problem.
2342
2343 :return: the right hand side of the problem
2344 :rtype: `Data`
2345 """
2346 return self.getSystem()[1]
2347
2348 def createRightHandSide(self):
2349 """
2350 Returns an instance of a new right hand side.
2351 """
2352 self.trace("New right hand side is allocated.")
2353 if self.getNumEquations()>1:
2354 return escore.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2355 else:
2356 return escore.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2357
2358 def createSolution(self):
2359 """
2360 Returns an instance of a new solution.
2361 """
2362 self.trace("New solution is allocated.")
2363 if self.getNumSolutions()>1:
2364 return escore.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2365 else:
2366 return escore.Data(0.,(),self.getFunctionSpaceForSolution(),True)
2367
2368 def resetSolution(self):
2369 """
2370 Sets the solution to zero.
2371 """
2372 if self.__solution.isEmpty():
2373 self.__solution=self.createSolution()
2374 else:
2375 self.__solution.setToZero()
2376 self.trace("Solution is reset to zero.")
2377
2378 def setSolution(self,u, validate=True):
2379 """
2380 Sets the solution assuming that makes the system valid with the tolrance
2381 defined by the solver options
2382 """
2383 if validate:
2384 self.__solution_rtol=self.getSolverOptions().getTolerance()
2385 self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2386 self.validSolution()
2387 self.__solution=u
2388 def getCurrentSolution(self):
2389 """
2390 Returns the solution in its current state.
2391 """
2392 if self.__solution.isEmpty(): self.__solution=self.createSolution()
2393 return self.__solution
2394
2395 def resetRightHandSide(self):
2396 """
2397 Sets the right hand side to zero.
2398 """
2399 if self.__righthandside.isEmpty():
2400 self.__righthandside=self.createRightHandSide()
2401 else:
2402 self.__righthandside.setToZero()
2403 self.trace("Right hand side is reset to zero.")
2404
2405 def getCurrentRightHandSide(self):
2406 """
2407 Returns the right hand side in its current state.
2408 """
2409 if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2410 return self.__righthandside
2411
2412 def resetOperator(self):
2413 """
2414 Makes sure that the operator is instantiated and returns it initialized
2415 with zeros.
2416 """
2417 if self.getOperatorType() == None:
2418 if self.isUsingLumping():
2419 self.__operator=self.createSolution()
2420 else:
2421 self.__operator=self.createOperator()
2422 self.__operator_type=self.getRequiredOperatorType()
2423 else:
2424 if self.isUsingLumping():
2425 self.__operator.setToZero()
2426 else:
2427 if self.getOperatorType() == self.getRequiredOperatorType():
2428 self.__operator.resetValues()
2429 else:
2430 self.__operator=self.createOperator()
2431 self.__operator_type=self.getRequiredOperatorType()
2432 self.trace("Operator reset to zero")
2433
2434 def getCurrentOperator(self):
2435 """
2436 Returns the operator in its current state.
2437 """
2438 return self.__operator
2439
2440 def setValue(self,**coefficients):
2441 """
2442 Sets new values to coefficients.
2443
2444 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2445 """
2446 # check if the coefficients are legal:
2447 for i in list(coefficients.keys()):
2448 if not self.hasCoefficient(i):
2449 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%i)
2450 # if the number of unknowns or equations is still unknown we try to estimate them:
2451 if self.__numEquations==None or self.__numSolutions==None:
2452 for i,d in list(coefficients.items()):
2453 if hasattr(d,"shape"):
2454 s=d.shape
2455 elif hasattr(d,"getShape"):
2456 s=d.getShape()
2457 else:
2458 s=numpy.array(d).shape
2459 if s!=None:
2460 # get number of equations and number of unknowns:
2461 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2462 if res==None:
2463 raise IllegalCoefficientValue("Illegal shape %s of coefficient %s"%(s,i))
2464 else:
2465 if self.__numEquations==None: self.__numEquations=res[0]
2466 if self.__numSolutions==None: self.__numSolutions=res[1]
2467 if self.__numEquations==None: raise UndefinedPDEError("unidentified number of equations")
2468 if self.__numSolutions==None: raise UndefinedPDEError("unidentified number of solutions")
2469 # now we check the shape of the coefficient if numEquations and numSolutions are set:
2470 for i,d in list(coefficients.items()):
2471 try:
2472 self.__COEFFICIENTS[i].setValue(self.getDomain(),
2473 self.getNumEquations(),self.getNumSolutions(),
2474 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2475 self.alteredCoefficient(i)
2476 except IllegalCoefficientFunctionSpace as m:
2477 # if the function space is wrong then we try the reduced version:
2478 i_red=i+"_reduced"
2479 if (not i_red in list(coefficients.keys())) and i_red in list(self.__COEFFICIENTS.keys()):
2480 try:
2481 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2482 self.getNumEquations(),self.getNumSolutions(),
2483 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2484 self.alteredCoefficient(i_red)
2485 except IllegalCoefficientValue as m:
2486 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2487 except IllegalCoefficientFunctionSpace as m:
2488 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2489 else:
2490 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2491 except IllegalCoefficientValue as m:
2492 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2493 self.__altered_coefficients=True
2494
2495 # ==========================================================================
2496 # methods that are typically overwritten when implementing a particular
2497 # linear problem
2498 # ==========================================================================
2499 def getRequiredOperatorType(self):
2500 """
2501 Returns the system type which needs to be used by the current set up.
2502
2503 :note: Typically this method is overwritten when implementing a
2504 particular linear problem.
2505 """
2506 return None
2507
2508 def createOperator(self):
2509 """
2510 Returns an instance of a new operator.
2511
2512 :note: This method is overwritten when implementing a particular
2513 linear problem.
2514 """
2515 return escore.Operator()
2516
2517 def checkSymmetry(self,verbose=True):
2518 """
2519 Tests the PDE for symmetry.
2520
2521 :param verbose: if set to True or not present a report on coefficients
2522 which break the symmetry is printed
2523 :type verbose: ``bool``
2524 :return: True if the problem is symmetric
2525 :rtype: ``bool``
2526 :note: Typically this method is overwritten when implementing a
2527 particular linear problem.
2528 """
2529 out=True
2530 return out
2531
2532 def getSolution(self,**options):
2533 """
2534 Returns the solution of the problem.
2535
2536 :return: the solution
2537 :rtype: `Data`
2538
2539 :note: This method is overwritten when implementing a particular
2540 linear problem.
2541 """
2542 return self.getCurrentSolution()
2543
2544 def getSystem(self):
2545 """
2546 Returns the operator and right hand side of the PDE.
2547
2548 :return: the discrete version of the PDE
2549 :rtype: ``tuple`` of `Operator` and `Data`.
2550
2551 :note: This method is overwritten when implementing a particular
2552 linear problem.
2553 """
2554 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2555
2556 class LinearPDE(LinearProblem):
2557 """
2558 This class is used to define a general linear, steady, second order PDE
2559 for an unknown function *u* on a given domain defined through a
2560 `Domain` object.
2561
2562 For a single PDE having a solution with a single component the linear PDE
2563 is defined in the following form:
2564
2565 *-(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)*
2566
2567 where *grad(F)* denotes the spatial derivative of *F*. Einstein's
2568 summation convention, ie. summation over indexes appearing twice in a term
2569 of a sum performed, is used.
2570 The coefficients *A*, *B*, *C*, *D*, *X* and *Y* have to be specified
2571 through `Data` objects in `Function` and
2572 the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
2573 *X_reduced* and *Y_reduced* have to be specified through
2574 `Data` objects in `ReducedFunction`.
2575 It is also allowed to use objects that can be converted into such
2576 `Data` objects. *A* and *A_reduced* are rank two, *B*,
2577 *C*, *X*, *B_reduced*, *C_reduced* and *X_reduced* are rank one and
2578 *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
2579
2580 The following natural boundary conditions are considered:
2581
2582 *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*
2583
2584 where *n* is the outer normal field. Notice that the coefficients *A*,
2585 *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
2586 PDE. The coefficients *d* and *y* are each a scalar in
2587 `FunctionOnBoundary` and the coefficients
2588 *d_reduced* and *y_reduced* are each a scalar in
2589 `ReducedFunctionOnBoundary`.
2590
2591 Constraints for the solution prescribe the value of the solution at certain
2592 locations in the domain. They have the form
2593
2594 *u=r* where *q>0*
2595
2596 *r* and *q* are each scalar where *q* is the characteristic function
2597 defining where the constraint is applied. The constraints override any
2598 other condition set by the PDE or the boundary condition.
2599
2600 The PDE is symmetrical if
2601
2602 *A[i,j]=A[j,i]* and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]*
2603 and *B_reduced[j]=C_reduced[j]*
2604
2605 For a system of PDEs and a solution with several components the PDE has the
2606 form
2607
2608 *-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]*
2609
2610 *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
2611 *C_reduced* are each of rank three, *D*, *D_reduced*, *X_reduced* and
2612 *X* are each of rank two and *Y* and *Y_reduced* are of rank one.
2613 The natural boundary conditions take the form:
2614
2615 *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]*
2616
2617 The coefficient *d* is of rank two and *y* is of rank one both in
2618 `FunctionOnBoundary`. The coefficients
2619 *d_reduced* is of rank two and *y_reduced* is of rank one both in
2620 `ReducedFunctionOnBoundary`.
2621
2622 Constraints take the form
2623
2624 *u[i]=r[i]* where *q[i]>0*
2625
2626 *r* and *q* are each rank one. Notice that at some locations not
2627 necessarily all components must have a constraint.
2628
2629 The system of PDEs is symmetrical if
2630
2631 - *A[i,j,k,l]=A[k,l,i,j]*
2632 - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
2633 - *B[i,j,k]=C[k,i,j]*
2634 - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2635 - *D[i,k]=D[i,k]*
2636 - *D_reduced[i,k]=D_reduced[i,k]*
2637 - *d[i,k]=d[k,i]*
2638 - *d_reduced[i,k]=d_reduced[k,i]*
2639
2640 `LinearPDE` also supports solution discontinuities over a contact region
2641 in the domain. To specify the conditions across the discontinuity we are
2642 using the generalised flux *J* which, in the case of a system of PDEs
2643 and several components of the solution, is defined as
2644
2645 *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]*
2646
2647 For the case of single solution component and single PDE *J* is defined as
2648
2649 *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]*
2650
2651 In the context of discontinuities *n* denotes the normal on the
2652 discontinuity pointing from side 0 towards side 1 calculated from
2653 `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2654 For a system of PDEs the contact condition takes the form
2655
2656 *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]*
2657
2658 where *J0* and *J1* are the fluxes on side 0 and side 1 of the
2659 discontinuity, respectively. *jump(u)*, which is the difference of the
2660 solution at side 1 and at side 0, denotes the jump of *u* across
2661 discontinuity along the normal calculated by `jump`.
2662 The coefficient *d_contact* is of rank two and *y_contact* is of rank one
2663 both in `FunctionOnContactZero` or
2664 `FunctionOnContactOne`.
2665 The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
2666 is of rank one both in `ReducedFunctionOnContactZero`
2667 or `ReducedFunctionOnContactOne`.
2668 In case of a single PDE and a single component solution the contact
2669 condition takes the form
2670
2671 *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
2672
2673 In this case the coefficient *d_contact* and *y_contact* are each scalar
2674 both in `FunctionOnContactZero` or
2675 `FunctionOnContactOne` and the coefficient
2676 *d_contact_reduced* and *y_contact_reduced* are each scalar both in
2677 `ReducedFunctionOnContactZero` or
2678 `ReducedFunctionOnContactOne`.
2679
2680 Typical usage::
2681
2682 p = LinearPDE(dom)
2683 p.setValue(A=kronecker(dom), D=1, Y=0.5)
2684 u = p.getSolution()
2685
2686 """
2687
2688 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2689 """
2690 Initializes a new linear PDE.
2691
2692 :param domain: domain of the PDE
2693 :type domain: `Domain`
2694 :param numEquations: number of equations. If ``None`` the number of
2695 equations is extracted from the PDE coefficients.
2696 :param numSolutions: number of solution components. If ``None`` the number
2697 of solution components is extracted from the PDE
2698 coefficients.
2699 :param debug: if True debug information is printed
2700
2701 """
2702 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2703 #
2704 # the coefficients of the PDE:
2705 #
2706 self.introduceCoefficients(
2707 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2708 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2709 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2710 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2711 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2712 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2713 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2714 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2715 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2716 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2717 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2718 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2719 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2720 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2721 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2722 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2723 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2724 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2725 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2726 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2727 d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2728 y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2729 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2730 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2731
2732 def __str__(self):
2733 """
2734 Returns the string representation of the PDE.
2735
2736 :return: a simple representation of the PDE
2737 :rtype: ``str``
2738 """
2739 return "<LinearPDE %d>"%id(self)
2740
2741 def getRequiredOperatorType(self):
2742 """
2743 Returns the system type which needs to be used by the current set up.
2744 """
2745 if self.isUsingLumping():
2746 return "__ESCRIPT_DATA"
2747 else:
2748 solver_options=self.getSolverOptions()
2749 return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2750
2751 def checkSymmetry(self,verbose=True):
2752 """
2753 Tests the PDE for symmetry.
2754
2755 :param verbose: if set to True or not present a report on coefficients
2756 which break the symmetry is printed.
2757 :type verbose: ``bool``
2758 :return: True if the PDE is symmetric
2759 :rtype: `bool`
2760 :note: This is a very expensive operation. It should be used for
2761 degugging only! The symmetry flag is not altered.
2762 """
2763 out=True
2764 out=out and self.checkSymmetricTensor("A", verbose)
2765 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2766 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2767 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2768 out=out and self.checkSymmetricTensor("D", verbose)
2769 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2770 out=out and self.checkSymmetricTensor("d", verbose)
2771 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2772 out=out and self.checkSymmetricTensor("d_contact", verbose)
2773 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2774 out=out and self.checkSymmetricTensor("d_dirac", verbose)
2775 return out
2776
2777 def createOperator(self):
2778 """
2779 Returns an instance of a new operator.
2780 """
2781 optype=self.getRequiredOperatorType()
2782 self.trace("New operator of type %s is allocated."%optype)
2783 return self.getDomain().newOperator( \
2784 self.getNumEquations(), \
2785 self.getFunctionSpaceForEquation(), \
2786 self.getNumSolutions(), \
2787 self.getFunctionSpaceForSolution(), \
2788 optype)
2789
2790 def getSolution(self):
2791 """
2792 Returns the solution of the PDE.
2793
2794 :return: the solution
2795 :rtype: `Data`
2796 """
2797 option_class=self.getSolverOptions()
2798 if not self.isSolutionValid():
2799 mat,f=self.getSystem()
2800 if self.isUsingLumping():
2801 if not util.inf(abs(mat)) > 0.:
2802 raise ZeroDivisionError("Lumped mass matrix has zero entry (try order 1 elements or HRZ lumping).")
2803 self.setSolution(f*1/mat)
2804 else:
2805 self.trace("PDE is resolved.")
2806 self.trace("solver options: %s"%str(option_class))
2807 self.setSolution(mat.solve(f,option_class))
2808 return self.getCurrentSolution()
2809
2810 def getSystem(self):
2811 """
2812 Returns the operator and right hand side of the PDE.
2813
2814 :return: the discrete version of the PDE
2815 :rtype: ``tuple`` of `Operator` and
2816 `Data`
2817 """
2818 if not self.isOperatorValid() or not self.isRightHandSideValid():
2819 if self.isUsingLumping():
2820 if not self.isOperatorValid():
2821 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2822 raise TypeError("Lumped matrix requires same order for equations and unknowns")
2823 if not self.getCoefficient("A").isEmpty():
2824 raise ValueError("coefficient A in lumped matrix may not be present.")
2825 if not self.getCoefficient("B").isEmpty():
2826 raise ValueError("coefficient B in lumped matrix may not be present.")
2827 if not self.getCoefficient("C").isEmpty():
2828 raise ValueError("coefficient C in lumped matrix may not be present.")
2829 if not self.getCoefficient("d_contact").isEmpty():
2830 raise ValueError("coefficient d_contact in lumped matrix may not be present.")
2831 if not self.getCoefficient("A_reduced").isEmpty():
2832 raise ValueError("coefficient A_reduced in lumped matrix may not be present.")
2833 if not self.getCoefficient("B_reduced").isEmpty():
2834 raise ValueError("coefficient B_reduced in lumped matrix may not be present.")
2835 if not self.getCoefficient("C_reduced").isEmpty():
2836 raise ValueError("coefficient C_reduced in lumped matrix may not be present.")
2837 if not self.getCoefficient("d_contact_reduced").isEmpty():
2838 raise ValueError("coefficient d_contact_reduced in lumped matrix may not be present.")
2839 D=self.getCoefficient("D")
2840 d=self.getCoefficient("d")
2841 D_reduced=self.getCoefficient("D_reduced")
2842 d_reduced=self.getCoefficient("d_reduced")
2843 d_dirac=self.getCoefficient("d_dirac")
2844
2845 if not D.isEmpty():
2846 if self.getNumSolutions()>1:
2847 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2848 else:
2849 D_times_e=D
2850 else:
2851 D_times_e=escore.Data()
2852 if not d.isEmpty():
2853 if self.getNumSolutions()>1:
2854 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2855 else:
2856 d_times_e=d
2857 else:
2858 d_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_reduced.isEmpty():
2869 if self.getNumSolutions()>1:
2870 d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2871 else:
2872 d_reduced_times_e=d_reduced
2873 else:
2874 d_reduced_times_e=escore.Data()
2875
2876 if not d_dirac.isEmpty():
2877 if self.getNumSolutions()>1:
2878 d_dirac_times_e=util.matrix_mult(d_dirac,numpy.ones((self.getNumSolutions(),)))
2879 else:
2880 d_reduced_dirac_e=d_dirac
2881 else:
2882 d_dirac_times_e=escore.Data()
2883
2884 self.resetOperator()
2885 operator=self.getCurrentOperator()
2886 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2887 hrz_lumping=( self.getSolverOptions().getSolverMethod() == SolverOptions.HRZ_LUMPING )
2888 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e, hrz_lumping )
2889 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, escore.Data(), hrz_lumping)
2890 else:
2891 self.getDomain().addPDEToRHS(operator, \
2892 escore.Data(), \
2893 D_times_e, \
2894 d_times_e,\
2895 escore.Data(),\
2896 d_dirac_times_e)
2897 self.getDomain().addPDEToRHS(operator, \
2898 escore.Data(), \
2899 D_reduced_times_e, \
2900 d_reduced_times_e,\
2901 escore.Data(), \
2902 escore.Data())
2903 self.trace("New lumped operator has been built.")
2904 if not self.isRightHandSideValid():
2905 self.resetRightHandSide()
2906 righthandside=self.getCurrentRightHandSide()
2907 self.getDomain().addPDEToRHS(righthandside, \
2908 self.getCoefficient("X"), \
2909 self.getCoefficient("Y"),\
2910 self.getCoefficient("y"),\
2911 self.getCoefficient("y_contact"), \
2912 self.getCoefficient("y_dirac"))
2913 self.getDomain().addPDEToRHS(righthandside, \
2914 self.getCoefficient("X_reduced"), \
2915 self.getCoefficient("Y_reduced"),\
2916 self.getCoefficient("y_reduced"),\
2917 self.getCoefficient("y_contact_reduced"), \
2918 escore.Data())
2919 self.trace("New right hand side has been built.")
2920 self.validRightHandSide()
2921 self.insertConstraint(rhs_only=False)
2922 self.validOperator()
2923 else:
2924 if not self.isOperatorValid() and not self.isRightHandSideValid():
2925 self.resetRightHandSide()
2926 righthandside=self.getCurrentRightHandSide()
2927 self.resetOperator()
2928 operator=self.getCurrentOperator()
2929 self.getDomain().addPDEToSystem(operator,righthandside, \
2930 self.getCoefficient("A"), \
2931 self.getCoefficient("B"), \
2932 self.getCoefficient("C"), \
2933 self.getCoefficient("D"), \
2934 self.getCoefficient("X"), \
2935 self.getCoefficient("Y"), \
2936 self.getCoefficient("d"), \
2937 self.getCoefficient("y"), \
2938 self.getCoefficient("d_contact"), \
2939 self.getCoefficient("y_contact"), \
2940 self.getCoefficient("d_dirac"), \
2941 self.getCoefficient("y_dirac"))
2942 self.getDomain().addPDEToSystem(operator,righthandside, \
2943 self.getCoefficient("A_reduced"), \
2944 self.getCoefficient("B_reduced"), \
2945 self.getCoefficient("C_reduced"), \
2946 self.getCoefficient("D_reduced"), \
2947 self.getCoefficient("X_reduced"), \
2948 self.getCoefficient("Y_reduced"), \
2949 self.getCoefficient("d_reduced"), \
2950 self.getCoefficient("y_reduced"), \
2951 self.getCoefficient("d_contact_reduced"), \
2952 self.getCoefficient("y_contact_reduced"), \
2953 escore.Data(), \
2954 escore.Data())
2955 self.insertConstraint(rhs_only=False)
2956 self.trace("New system has been built.")
2957 self.validOperator()
2958 self.validRightHandSide()
2959 elif not self.isRightHandSideValid():
2960 self.resetRightHandSide()
2961 righthandside=self.getCurrentRightHandSide()
2962 self.getDomain().addPDEToRHS(righthandside,
2963 self.getCoefficient("X"), \
2964 self.getCoefficient("Y"),\
2965 self.getCoefficient("y"),\
2966 self.getCoefficient("y_contact"), \
2967 self.getCoefficient("y_dirac") )
2968 self.getDomain().addPDEToRHS(righthandside,
2969 self.getCoefficient("X_reduced"), \
2970 self.getCoefficient("Y_reduced"),\
2971 self.getCoefficient("y_reduced"),\
2972 self.getCoefficient("y_contact_reduced"), \
2973 escore.Data())
2974 self.insertConstraint(rhs_only=True)
2975 self.trace("New right hand side has been built.")
2976 self.validRightHandSide()
2977 elif not self.isOperatorValid():
2978 self.resetOperator()
2979 operator=self.getCurrentOperator()
2980 self.getDomain().addPDEToSystem(operator,escore.Data(), \
2981 self.getCoefficient("A"), \
2982 self.getCoefficient("B"), \
2983 self.getCoefficient("C"), \
2984 self.getCoefficient("D"), \
2985 escore.Data(), \
2986 escore.Data(), \
2987 self.getCoefficient("d"), \
2988 escore.Data(),\
2989 self.getCoefficient("d_contact"), \
2990 escore.Data(), \
2991 self.getCoefficient("d_dirac"), \
2992 escore.Data())
2993 self.getDomain().addPDEToSystem(operator,escore.Data(), \
2994 self.getCoefficient("A_reduced"), \
2995 self.getCoefficient("B_reduced"), \
2996 self.getCoefficient("C_reduced"), \
2997 self.getCoefficient("D_reduced"), \
2998 escore.Data(), \
2999 escore.Data(), \
3000 self.getCoefficient("d_reduced"), \
3001 escore.Data(),\
3002 self.getCoefficient("d_contact_reduced"), \
3003 escore.Data(), \
3004 escore.Data(), \
3005 escore.Data())
3006 self.insertConstraint(rhs_only=False)
3007 self.trace("New operator has been built.")
3008 self.validOperator()
3009 self.setSystemStatus()
3010 self.trace("System status is %s."%self.getSystemStatus())
3011 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
3012
3013 def insertConstraint(self, rhs_only=False):
3014 """
3015 Applies the constraints defined by q and r to the PDE.
3016
3017 :param rhs_only: if True only the right hand side is altered by the
3018 constraint
3019 :type rhs_only: ``bool``
3020 """
3021 q=self.getCoefficient("q")
3022 r=self.getCoefficient("r")
3023 righthandside=self.getCurrentRightHandSide()
3024 operator=self.getCurrentOperator()
3025
3026 if not q.isEmpty():
3027 if r.isEmpty():
3028 r_s=self.createSolution()
3029 else:
3030 r_s=r
3031 if not rhs_only and not operator.isEmpty():
3032 if self.isUsingLumping():
3033 operator.copyWithMask(escore.Data(1.,q.getShape(),q.getFunctionSpace()),q)
3034 else:
3035 row_q=escore.Data(q,self.getFunctionSpaceForEquation())
3036 col_q=escore.Data(q,self.getFunctionSpaceForSolution())
3037 u=self.createSolution()
3038 u.copyWithMask(r_s,col_q)
3039 righthandside-=operator*u
3040 operator.nullifyRowsAndCols(row_q,col_q,1.)
3041 righthandside.copyWithMask(r_s,q)
3042
3043 def setValue(self,**coefficients):
3044 """
3045 Sets new values to coefficients.
3046
3047 :param coefficients: new values assigned to coefficients
3048 :keyword A: value for coefficient ``A``
3049 :type A: any type that can be cast to a `Data` object on
3050 `Function`
3051 :keyword A_reduced: value for coefficient ``A_reduced``
3052 :type A_reduced: any type that can be cast to a `Data`
3053 object on `ReducedFunction`
3054 :keyword B: value for coefficient ``B``
3055 :type B: any type that can be cast to a `Data` object on
3056 `Function`
3057 :keyword B_reduced: value for coefficient ``B_reduced``
3058 :type B_reduced: any type that can be cast to a `Data`
3059 object on `ReducedFunction`
3060 :keyword C: value for coefficient ``C``
3061 :type C: any type that can be cast to a `Data` object on
3062 `Function`
3063 :keyword C_reduced: value for coefficient ``C_reduced``
3064 :type C_reduced: any type that can be cast to a `Data`
3065 object on `ReducedFunction`
3066 :keyword D: value for coefficient ``D``
3067 :type D: any type that can be cast to a `Data` object on
3068 `Function`
3069 :keyword D_reduced: value for coefficient ``D_reduced``
3070 :type D_reduced: any type that can be cast to a `Data`
3071 object on `ReducedFunction`
3072 :keyword X: value for coefficient ``X``
3073 :type X: any type that can be cast to a `Data` object on
3074 `Function`
3075 :keyword X_reduced: value for coefficient ``X_reduced``
3076 :type X_reduced: any type that can be cast to a `Data`
3077 object on `ReducedFunction`
3078 :keyword Y: value for coefficient ``Y``
3079 :type Y: any type that can be cast to a `Data` object on
3080 `Function`
3081 :keyword Y_reduced: value for coefficient ``Y_reduced``
3082 :type Y_reduced: any type that can be cast to a `Data`
3083 object on `ReducedFunction`
3084 :keyword d: value for coefficient ``d``
3085 :type d: any type that can be cast to a `Data` object on
3086 `FunctionOnBoundary`
3087 :keyword d_reduced: value for coefficient ``d_reduced``
3088 :type d_reduced: any type that can be cast to a `Data`
3089 object on `ReducedFunctionOnBoundary`
3090 :keyword y: value for coefficient ``y``
3091 :type y: any type that can be cast to a `Data` object on
3092 `FunctionOnBoundary`
3093 :keyword d_contact: value for coefficient ``d_contact``
3094 :type d_contact: any type that can be cast to a `Data`
3095 object on `FunctionOnContactOne`
3096 or `FunctionOnContactZero`
3097 :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
3098 :type d_contact_reduced: any type that can be cast to a `Data`
3099 object on `ReducedFunctionOnContactOne`
3100 or `ReducedFunctionOnContactZero`
3101 :keyword y_contact: value for coefficient ``y_contact``
3102 :type y_contact: any type that can be cast to a `Data`
3103 object on `FunctionOnContactOne`
3104 or `FunctionOnContactZero`
3105 :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
3106 :type y_contact_reduced: any type that can be cast to a `Data`
3107 object on `ReducedFunctionOnContactOne`
3108 or `ReducedFunctionOnContactZero`
3109 :keyword d_dirac: value for coefficient ``d_dirac``
3110 :type d_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3111 :keyword y_dirac: value for coefficient ``y_dirac``
3112 :type y_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
3113 :keyword r: values prescribed to the solution at the locations of
3114 constraints
3115 :type r: any type that can be cast to a `Data` object on
3116 `Solution` or `ReducedSolution`
3117 depending on whether reduced order is used for the solution
3118 :keyword q: mask for location of constraints
3119 :type q: any type that can be cast to a `Data` object on
3120 `Solution` or `ReducedSolution`
3121 depending on whether reduced order is used for the
3122 representation of the equation
3123 :raise IllegalCoefficient: if an unknown coefficient keyword is used
3124 """
3125 super(LinearPDE,self).setValue(**coefficients)
3126 # check if the systrem is inhomogeneous:
3127 if len(coefficients)>0 and not self.isUsingLumping():
3128 q=self.getCoefficient("q")
3129 r=self.getCoefficient("r")
3130 if not q.isEmpty() and not r.isEmpty():
3131 if util.Lsup(q*r)>0.:
3132 self.trace("Inhomogeneous constraint detected.")
3133 self.invalidateSystem()
3134
3135
3136 def getResidual(self,u=None):
3137 """
3138 Returns the residual of u or the current solution if u is not present.
3139
3140 :param u: argument in the residual calculation. It must be representable
3141 in `self.getFunctionSpaceForSolution()`. If u is not present
3142 or equals ``None`` the current solution is used.
3143 :type u: `Data` or None
3144 :return: residual of u
3145 :rtype: `Data`
3146 """
3147 if u==None:
3148 return self.getOperator()*self.getSolution()-self.getRightHandSide()
3149 else:
3150 return self.getOperator()*escore.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
3151
3152 def getFlux(self,u=None):
3153 """
3154 Returns the flux *J* for a given *u*.
3155
3156 *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]*
3157
3158 or
3159
3160 *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]*
3161
3162 :param u: argument in the flux. If u is not present or equals `None` the
3163 current solution is used.
3164 :type u: `Data` or None
3165 :return: flux
3166 :rtype: `Data`
3167 """
3168 if u==None: u=self.getSolution()
3169 if self.getNumEquations()>1:
3170 out = escore.Data(0.,(self.getNumEquations(),self.getDim()),self.getFunctionSpaceForCoefficient("X"))
3171 else:
3172 out = escore.Data(0.,(self.getDim(), ),self.getFunctionSpaceForCoefficient("X"))
3173
3174 A=self.getCoefficient("A")
3175 if not A.isEmpty():
3176 out+=util.tensormult(A,util.grad(u,self.getFunctionSpaceForCoefficient("A")))
3177
3178 B=self.getCoefficient("B")
3179 if not B.isEmpty():
3180 if B.getRank() == 1:
3181 out+=B * u
3182 else:
3183 out+=util.generalTensorProduct(B,u,axis_offset=1)
3184
3185 X=self.getCoefficient("X")
3186 if not X.isEmpty():
3187 out-=X
3188
3189 A_reduced=self.getCoefficient("A_reduced")
3190 if not A_reduced.isEmpty():
3191 out+=util.tensormult(A_reduced, util.grad(u,self.getFunctionSpaceForCoefficient("A_reduced"))) \
3192
3193 B_reduced=self.getCoefficient("B_reduced")
3194 if not B_reduced.isEmpty():
3195 if B_reduced.getRank() == 1:
3196 out+=B_reduced*u
3197 else:
3198 out+=util.generalTensorProduct(B_reduced,u,axis_offset=1)
3199
3200 X_reduced=self.getCoefficient("X_reduced")
3201 if not X_reduced.isEmpty():
3202 out-=X_reduced
3203 return out
3204
3205 class Poisson(LinearPDE):
3206 """
3207 Class to define a Poisson equation problem. This is generally a
3208 `LinearPDE` of the form
3209
3210 *-grad(grad(u)[j])[j] = f*
3211
3212 with natural boundary conditions
3213
3214 *n[j]*grad(u)[j] = 0*
3215
3216 and constraints:
3217
3218 *u=0* where *q>0*
3219
3220 """
3221
3222 def __init__(self,domain,debug=False):
3223 """
3224 Initializes a new Poisson equation.
3225
3226 :param domain: domain of the PDE
3227 :type domain: `Domain`
3228 :param debug: if True debug information is printed
3229
3230 """
3231 super(Poisson, self).__init__(domain,1,1,debug)
3232 self.introduceCoefficients(
3233 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3234 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
3235 self.setSymmetryOn()
3236
3237 def setValue(self,**coefficients):
3238 """
3239 Sets new values to coefficients.
3240
3241 :param coefficients: new values assigned to coefficients
3242 :keyword f: value for right hand side *f*
3243 :type f: any type that can be cast to a `Scalar` object
3244 on `Function`
3245 :keyword q: mask for location of constraints
3246 :type q: any type that can be cast to a rank zero `Data`
3247 object on `Solution` or
3248 `ReducedSolution` depending on whether
3249 reduced order is used for the representation of the equation
3250 :raise IllegalCoefficient: if an unknown coefficient keyword is used
3251 """
3252 super(Poisson, self).setValue(**coefficients)
3253
3254
3255 def getCoefficient(self,name):
3256 """
3257 Returns the value of the coefficient ``name`` of the general PDE.
3258
3259 :param name: name of the coefficient requested
3260 :type name: ``string``
3261 :return: the value of the coefficient ``name``
3262 :rtype: `Data`
3263 :raise IllegalCoefficient: invalid coefficient name
3264 :note: This method is called by the assembling routine to map the Poisson
3265 equation onto the general PDE.
3266 """
3267 if name == "A" :
3268 return escore.Data(util.kronecker(self.getDim()),escore.Function(self.getDomain()))
3269 elif name == "Y" :
3270 return self.getCoefficient("f")
3271 elif name == "Y_reduced" :
3272 return self.getCoefficient("f_reduced")
3273 else:
3274 return super(Poisson, self).getCoefficient(name)
3275
3276 class Helmholtz(LinearPDE):
3277 """
3278 Class to define a Helmholtz equation problem. This is generally a
3279 `LinearPDE` of the form
3280
3281 *omega*u - grad(k*grad(u)[j])[j] = f*
3282
3283 with natural boundary conditions
3284
3285 *k*n[j]*grad(u)[j] = g- alphau*
3286
3287 and constraints:
3288
3289 *u=r* where *q>0*
3290
3291 """
3292
3293 def __init__(self,domain,debug=False):
3294 """
3295 Initializes a new Helmholtz equation.
3296
3297 :param domain: domain of the PDE
3298 :type domain: `Domain`
3299 :param debug: if True debug information is printed
3300
3301 """
3302 super(Helmholtz, self).__init__(domain,1,1,debug)
3303 self.introduceCoefficients(
3304 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3305 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3306 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3307 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3308 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
3309 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
3310 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
3311 self.setSymmetryOn()
3312
3313 def setValue(self,**coefficients):
3314 """
3315 Sets new values to coefficients.
3316
3317 :param coefficients: new values assigned to coefficients
3318 :keyword omega: value for coefficient *omega*
3319 :type omega: any type that can be cast to a `Scalar`
3320 object on `Function`
3321 :keyword k: value for coefficient *k*
3322 :type k: any type that can be cast to a `Scalar` object
3323 on `Function`
3324 :keyword f: value for right hand side *f*
3325 :type f: any type that can be cast to a `Scalar` object
3326 on `Function`
3327 :keyword alpha: value for right hand side *alpha*
3328 :type alpha: any type that can be cast to a `Scalar`
3329 object on `FunctionOnBoundary`
3330 :keyword g: value for right hand side *g*
3331 :type g: any type that can be cast to a `Scalar` object
3332 on `FunctionOnBoundary`
3333 :keyword r: prescribed values *r* for the solution in constraints
3334 :type r: any type that can be cast to a `Scalar` object
3335 on `Solution` or
3336 `ReducedSolution` depending on whether
3337 reduced order is used for the representation of the equation
3338 :keyword q: mask for the location of constraints
3339 :type q: any type that can be cast to a `Scalar` object
3340 on `Solution` or
3341 `ReducedSolution` depending on whether
3342 reduced order is used for the representation of the equation
3343 :raise IllegalCoefficient: if an unknown coefficient keyword is used
3344 """
3345 super(Helmholtz, self).setValue(**coefficients)
3346
3347 def getCoefficient(self,name):
3348 """
3349 Returns the value of the coefficient ``name`` of the general PDE.
3350
3351 :param name: name of the coefficient requested
3352 :type name: ``string``
3353 :return: the value of the coefficient ``name``
3354 :rtype: `Data`
3355 :raise IllegalCoefficient: invalid name
3356 """
3357 if name == "A" :
3358 if self.getCoefficient("k").isEmpty():
3359 return escore.Data(numpy.identity(self.getDim()),escore.Function(self.getDomain()))
3360 else:
3361 return escore.Data(numpy.identity(self.getDim()),escore.Function(self.getDomain()))*self.getCoefficient("k")
3362 elif name == "D" :
3363 return self.getCoefficient("omega")
3364 elif name == "Y" :
3365 return self.getCoefficient("f")
3366 elif name == "d" :
3367 return self.getCoefficient("alpha")
3368 elif name == "y" :
3369 return self.getCoefficient("g")
3370 elif name == "Y_reduced" :
3371 return self.getCoefficient("f_reduced")
3372 elif name == "y_reduced" :
3373 return self.getCoefficient("g_reduced")
3374 else:
3375 return super(Helmholtz, self).getCoefficient(name)
3376
3377 class VTIWavePDE(LinearPDE):
3378 """
3379 A class specifically for waves, passes along values to native implementation
3380 to save computational time.
3381 """
3382 def __init__(self,domain,c,numEquations=None,numSolutions=None,debug=False):
3383 """
3384 Initializes a new linear PDE.
3385
3386 :param domain: domain of the PDE
3387 :type domain: `Domain`
3388 :param numEquations: number of equations. If ``None`` the number of
3389 equations is extracted from the PDE coefficients.
3390 :param numSolutions: number of solution components. If ``None`` the number
3391 of solution components is extracted from the PDE
3392 coefficients.
3393 :param debug: if True debug information is printed
3394
3395 """
3396 super(VTIWavePDE, self).__init__(domain,numEquations,numSolutions,debug)
3397 #
3398 # the coefficients of the PDE:
3399 #
3400 self.introduceCoefficients(
3401 A