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