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