Source code for pygimli.solver.solver

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Finite-element solver and utility functions."""
from copy import deepcopy

import numpy as np
import pygimli as pg


[docs] def parseDictKey_(key, markers): return parseMarkersDictKey(key, markers)
[docs] def parseMarkersDictKey(key, markers): """ Parse dictionary key of type str to marker list. Utility function to parse a dictionary key string into a valid list of markers containing in a given markers list. Parameters ---------- key: str | int Supported are - int: single markers - '*': all markers - 'm1': Single marker - 'm1,m2': Comma separated list - ':': Slice wildcard - 'start:stop:step': Slice like syntax markers: [int] List of integers, e.g., cell or boundary markers Returns ------- mas: [int] List of integers described by key """ markers = pg.unique(markers) mas = None if isinstance(key, str): if key == '*': return markers if ',' in key: mas = [int(k) for k in key.split(',')] elif ':' in key: sse = key.split(':') start = markers[0] stop = markers[-1] + 1 step = 1 if len(sse) > 0: try: start = int(sse[0]) except BaseException: pass if len(sse) > 1: try: stop = int(sse[1]) except BaseException: pass if len(sse) > 2: try: step = int(sse[2]) except BaseException: pass mas = list(range(start, stop, step)) else: mas = [int(key)] else: mas = [int(key)] return [m for m in mas if m in markers]
[docs] def boundaryIdsFromDictKey(mesh, key, outside=True): """Find all boundaries matching a dictionary key. Attributes ---------- mesh: :gimliapi:`GIMLI::Mesh` key: str|int Representation for boundary marker. Will be parsed by :py:mod:`pygimli.solver.solver.parseMarkersDictKey` outside: bool [True] Only select outside boundaries. Returns ------- dict: {marker, [boundary.id()]} """ mas = pg.solver.parseMarkersDictKey(key, mesh.boundaryMarkers()) ret = dict() for m in mas: for i in pg.find(mesh.boundaryMarkers() == m): if m not in ret: ret[m] = [] if outside is True and not mesh.boundary(i).outside(): continue ret[m].append(i) return ret
[docs] def cellValues(mesh, arg, **kwargs): """Get a value for each cell. Returns a array or vector of length mesh.cellCount() based on arg. The preferable arg is a dictionary for the cell marker and the appropriate cell value. The designated value can be calculated using a callable(cell, **kwargs), which is called on demand. Attributes ---------- mesh: :gimliapi:`GIMLI::Mesh` Used if arg is callable arg: float | int | complex | ndarray | iterable | callable | dict Argument to be parsed as cell data. If arg is a dictionary, its key will be interpreted as cell marker: Dictionary is key: value. Value can be float, int, complex or ndarray. The last for anistropic or elastic tensors. Key can be integer for cell marker or str, which will be interpreted as splice or list. See examples or `py:mod:pygimli.solver.parseMarkersDictKey`. Iterable of length mesh.nodeCount() to be interpolated to cell centers. userData: class Used if arg is callable Returns ------- ret: :gimliapi:`GIMLI::RVector` | ndarray(mesh.cellCount(), xx ) Array of desired length filled with the appropriate values. Examples -------- >>> import pygimli as pg >>> mesh = pg.createGrid(x=range(5)) >>> mesh.setCellMarkers([1, 1, 2, 2]) >>> print(mesh.cellCount()) 4 >>> print(pg.solver.cellValues(mesh, [1, 2, 3, 4])) [1, 2, 3, 4] >>> print(pg.solver.cellValues(mesh, {1:1.0, 2:10})) [1.0, 1.0, 10, 10] >>> print(pg.solver.cellValues(mesh, {':':2.0})) [2.0, 2.0, 2.0, 2.0] >>> print(pg.solver.cellValues(mesh, {'0:2':3.0})) [3.0, 3.0, None, None] >>> print(pg.solver.cellValues(mesh, np.ones(mesh.nodeCount()))) 4 [1.0, 1.0, 1.0, 1.0] >>> print(np.array(pg.solver.cellValues(mesh, {'1:3' : np.diag([1.0, 2.0])}))) [[[1. 0.] [0. 2.]] <BLANKLINE> [[1. 0.] [0. 2.]] <BLANKLINE> [[1. 0.] [0. 2.]] <BLANKLINE> [[1. 0.] [0. 2.]]] >>> print(np.array(pg.solver.cellValues(mesh, {':' : pg.core.CMatrix(2, 2)}))) [[[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] <BLANKLINE> [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] <BLANKLINE> [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]] <BLANKLINE> [[0.+0.j 0.+0.j] [0.+0.j 0.+0.j]]] >>> print(pg.solver.cellValues(mesh, {'1,2':1 + 1j*2.0})) [(1+2j), (1+2j), (1+2j), (1+2j)] >>> def cellVal(c, b=1): ... return c.center()[0]*b >>> t = pg.solver.cellValues(mesh, {':' : cellVal}) >>> print([t[c.id()](c) for c in mesh.cells()]) [0.5, 1.5, 2.5, 3.5] """ if isinstance(arg, dict): try: val = list(arg.values())[0] except BaseException: pg.error("Can't interpret empty dictionary:", arg) val = 1.0 ret = [None] * mesh.cellCount() for key, val in arg.items(): if isinstance(key, str): mas = parseDictKey_(key, mesh.cellMarkers()) for m in mas: for i in pg.find(mesh.cellMarkers() == m): ret[i] = val else: for i in pg.find(mesh.cellMarkers() == key): ret[i] = val return ret # if arg have already the correct size if hasattr(arg, '__len__'): if len(arg) == mesh.cellCount(): return arg if len(arg) == mesh.nodeCount(): return pg.interpolate(mesh, arg, mesh.cellCenters()) # if arg if scalar or global data type, ndarray or Matrix but not the right # size assume global tensor if isinstance(arg, np.ndarray) or \ isinstance(arg, pg.core.RMatrix) or \ isinstance(arg, pg.core.CMatrix) or \ isinstance(arg, float) or \ isinstance(arg, int) or \ isinstance(arg, complex): return [arg]*mesh.cellCount() return parseArgToArray(arg, nDof=mesh.cellCount(), mesh=mesh, **kwargs)
[docs] def parseArgToArray(arg, nDof, mesh=None, userData={}): """ Parse array related arguments to create a valid value array. Parameters ---------- arg : float | int | iterable | callable The target array value that will be converted to an array. If arg is a callable with it must fulfill: :: arg(cell|node|boundary, userData={}) Where MeshEntity is one of :gimliapi:`GIMLI::Cell` , :gimliapi:`GIMLI::Node` or :gimliapi:`GIMLI::Boundary` depending on nDof, where nDof is mesh.cellCount(), mesh.nodeCount() or mesh.boundaryCount(), respectively. nDof : int | [int] Desired array size. mesh : :gimliapi:`GIMLI::Mesh` Used if arg is callable userData : class Used if arg is callable Returns ------- ret : :gimliapi:`GIMLI::RVector` Array of desired length filled with the appropriate values. """ # pg.warn('check if obsolete: parseArgToArray') if not hasattr(nDof, '__len__'): nDof = [nDof] try: return pg.Vector(nDof[0], float(arg)) except BaseException: pass if hasattr(arg, '__len__'): if isinstance(arg, np.ndarray): if len(arg) == nDof[0]: return arg else: raise Exception('Given array does not have requested (' + str(nDof) + ') size (' + str(len(arg)) + ')') for n in nDof: if len(arg) == n: return arg try: # [marker, val] || [[marker, val]] return parseMapToCellArray(arg, mesh) except BaseException: raise Exception("Array 'arg' has the wrong size: " + str(len(arg)) + " != " + str(nDof)) elif hasattr(arg, '__call__'): ret = pg.Vector(nDof[0], 0.0) if not mesh: raise Exception("Please provide a mesh for the callable" "argument to parse ") if nDof[0] == mesh.nodeCount(): for n in mesh.nodes(): if userData: ret[n.id()] = arg(node=n, userData=userData) else: ret[n.id()] = arg(node=n) elif nDof[0] == mesh.cellCount(): for c in mesh.cells(): if userData: ret[c.id()] = arg(cell=c, userData=userData) else: ret[c.id()] = arg(cell=c) elif nDof[0] == mesh.boundaryCount(): for b in mesh.boundaries(): if userData: ret[b.id()] = arg(boundary=b, userData=userData) else: ret[b.id()] = arg(boundary=b) else: raise Exception("Cannot parse callable argument " + str(nDof) + " nodes: " + str(mesh.nodeCount()) + " cells: " + str(mesh.cellCount())) return ret raise Exception("Cannot parse argument type " + str(type(arg)))
[docs] def generateBoundaryValue(boundary, arg, time=0.0, userData={}, expectList=False, nCoeff=1): """Generate a value for the given Boundary. TODO ---- * support for complex vals Parameters ---------- boundary: :gimliapi:`GIMLI::Boundary` or list of .. The related boundary. expectList: bool[False] Allow list values for Robin BC. arg: convertible | iterable | callable or list of .. - convertible into float - iterable of minimum length = boundary.id() - callable generator function If arg is a callable it must fulfill: :: arg(boundary=:gimliapi:`GIMLI::Boundary`, time=0.0, userData={}) The callable function arg have to return appropriate values for all nodes of the boundary or one value for all nodes (scalar field only). Value can be scalar or vector field value, e.g., return force values for all nodes at a boundary to return an ndarray((nodes, dims)), e.g. 'lambda _b: np.array([[forc_x, forc_y, forc_z] for n in _b.nodes()]).T' Returns ------- val: [float] Value for all nodes of the boundary. """ val = 0. if callable(arg): kwargs = dict() if time != 0.0 and time is not None: kwargs['time'] = time if userData is not None and userData.keys(): kwargs['userData'] = userData try: # val(boundary, time=0, userData={}) val = arg(boundary, **kwargs) except BaseException as e: print(arg, "(", kwargs, ")") pg.critical("Wrong arguments for callback function.", e) elif hasattr(arg, '__len__'): if callable(arg[0]): kwargs = arg[1] if time != 0.0 and time is not None: kwargs['time'] = time val = arg[0](boundary=boundary, **kwargs) else: val = arg else: try: val = float(arg) except ValueError: print(arg, val) pg.error("can't create boundary values.") # transform val into list of length nodeCount if expectList is True: if np.array(val).ndim != 2: val = np.atleast_1d(val) if isinstance(boundary, pg.core.Node): return val if nCoeff == 1 and expectList is False: if isinstance(val, float): val = np.ones(boundary.nodeCount(), dtype=float) * val if len(val) != boundary.nodeCount(): print(val) pg.critical("Boundary value cannot be generated for nCoeff=1 val:", val) else: val = np.atleast_2d(val) # pg._y(val) if len(val) != boundary.nodeCount() or val.shape[1] != nCoeff: val = np.tile(val, (boundary.nodeCount(), 1)) return val
[docs] def parseArgPairToBoundaryArray(pair, mesh): """ Parse boundary related pair argument to create a list of [ :gimliapi:`GIMLI::Boundary`, value|callable ]. Parameters ---------- pair: tuple - [marker, arg] - [marker, [callable, *kwargs]] - [marker, [arg_x, arg_y, arg_z]] - [boundary, arg] - ['*', arg] - [node, arg] - [[marker, ...], arg] (REMOVE ME because of bad design) - [[boundary,...], arg] (REMOVE ME because of bad design) - [marker, callable, *kwargs] (REMOVE ME because of bad design) - [[marker, ...], callable, *kwargs] (REMOVE ME because of bad design) arg will be parsed by :py:mod:`pygimli.solver.solver.generateBoundaryValue` and distributed to each boundary. Callable functions will be executed at run time. '*' is interpreted as all boundary elements with one neighboring cell mesh: :gimliapi:`GIMLI::Mesh` Used to find boundaries by marker. Returns ------- bc: list() [:gimliapi:`GIMLI::Boundary`, value|callable] """ bc = [] bounds = [] if isinstance(pair[1], list): # [marker, [callable, *kwargs]] if callable(pair[1][0]): pair = [pair[0]] + pair[1] if pair[0] == '*': mesh.createNeighborInfos() for b in mesh.boundaries(): if b.leftCell() is not None and b.rightCell() is None: bounds.append(b) elif isinstance(pair[0], int): bounds = mesh.findBoundaryByMarker(pair[0]) elif isinstance(pair[0], pg.core.Node): bc.append(pair) return bc for b in bounds: val = None if len(pair) > 2: val = pair[1:] else: val = pair[1] bc.append([b, val]) # print('-'*50) # print(b, pair[1], callable(pair[1])) # print('+'*50) # if callable(pair[1]): # # don't execute the callable here # # we want to call them at runtime # if len(pair) > 2: # val = pair[1:] # else: # val = pair[1] # else: # this will be executed # val = generateBoundaryValue(b, pair[1]) # print('#'*30) return bc
[docs] def parseArgToBoundaries(args, mesh): """ Parse boundary related arguments to create a valid boundary value list: [ :gimliapi:`GIMLI::Boundary`, value|callable ] TODO ---- - callable dynamic at runtime Parameters ---------- args : dict, float, callable Dictionary is preferred (key=value|callable). If args is just a callable or float every outer boundary is processed with args. List pairs will be removed or not correct parsed for vector valued problems. Callable will be evaluated at runtime. See examples. Else see :py:mod:`pygimli.solver.solver.parseArgPairToBoundaryArray` mesh : :gimliapi:`GIMLI::Mesh` Used to find boundaries by marker Returns ------- boundaries : list() [ :gimliapi:`GIMLI::Boundary`, value|callable ] Examples -------- >>> # no need to import matplotlib. pygimli show does >>> import pygimli as pg >>> import pygimli.meshtools as mt >>> plc = mt.createWorld([0, 0], [1, -1], worldMarker=0) >>> ax, _ = pg.show(plc, boundaryMarker=True) >>> mesh = mt.createMesh(plc) >>> # all four outer boundaries get value = 1.0 >>> b = pg.solver.parseArgToBoundaries(1.0, mesh) >>> print(len(b)) 4 >>> # all edges with marker 1 get value = 1.0 >>> b = pg.solver.parseArgToBoundaries({1: 1.0}, mesh) >>> print(len(b)) 1 >>> # same as above with marker 2 get value 2 >>> b = pg.solver.parseArgToBoundaries({'1': 1.0, 2 : 2.0}, mesh) >>> print(len(b)) 2 >>> # same as above with marker 3 get value 3 >>> b = pg.solver.parseArgToBoundaries({1:1., 2:2., 3:3.}, mesh) >>> print(len(b)) 3 >>> # Boundary values for vector valued problem >>> b = pg.solver.parseArgToBoundaries({1:[1.0, 1.0]}, mesh) >>> print(len(b), b[0][1]) 1 [1.0, 1.0] >>> # edges with marker 1 and 2 get value 1 >>> b = pg.solver.parseArgToBoundaries({'1,2':1.0}, mesh) >>> print(len(b)) 2 >>> b = pg.solver.parseArgToBoundaries({'1, 2, 3': 1.0}, mesh) >>> print(len(b)) 3 >>> b = pg.solver.parseArgToBoundaries({'1:4':1.0, 4:4.0}, mesh) >>> print(len(b)) 4 >>> b = pg.solver.parseArgToBoundaries({mesh.node(0):0.0}, mesh) >>> print(len(b)) 1 >>> def bCall(boundary): ... u = [] ... for i, n in enumerate(boundary.nodes()): ... u.append(i) ... return u >>> b = pg.solver.parseArgToBoundaries({1:bCall}, mesh) >>> print(len(b),b[0][1](b[0][0])) 1 [0, 1] >>> def bCall(boundary, a1, a2): ... return a1 + a2 >>> b = pg.solver.parseArgToBoundaries({1: [bCall, {'a1':2, 'a2':3}]}, mesh) >>> print(len(b), b[0][1][0](b[0][0], **b[0][1][1])) 1 5 >>> b = pg.solver.parseArgToBoundaries({1: [bCall, {'a1':1, 'a2':2}], ... 2: [bCall, {'a1':3, 'a2':4}]}, mesh) >>> print(len(b), b[0][1][0](b[0][0], **b[0][1][1])) 2 3 >>> b = pg.solver.parseArgToBoundaries({'1,2': [bCall, {'a1':4, 'a2':5}]}, mesh) >>> print(len(b), b[1][1][0](b[1][0], **b[1][1][1])) 2 9 >>> pg.wait() """ boundaries = list() if isinstance(args, dict): # try: # val = list(args.values())[0] # except BaseException as _: # return boundaries # pg.error("Can't interpret empty dictionary:", args) for key, val in args.items(): if isinstance(key, pg.core.stdVectorNodes) or \ isinstance(key, list) and isinstance(key[0], pg.core.Node): for n in key: boundaries += parseArgPairToBoundaryArray([n, val], mesh) elif isinstance(key, str) and key != '*': markers = parseDictKey_(key, mesh.boundaryMarkers()) for m in markers: boundaries += parseArgPairToBoundaryArray([m, val], mesh) else: boundaries += parseArgPairToBoundaryArray([key, val], mesh) return boundaries if hasattr(args, '__call__') or isinstance(args, float) or \ isinstance(args, int): return parseArgToBoundaries({'*': args}, mesh) else: raise Exception('cannot interpret boundary token', args) return boundaries
def _bcIsForVectorValues(bc, mesh): """Guess if boundary condition is supposed to be for vector valued problems """ verbose = False def testForV3(t): if verbose: print("test for v3", t) try: if callable(t): test = t(mesh.boundary(0)) else: test = t if verbose: print("test for v3 test", test) if hasattr(test, '__iter__'): test = np.array(test) # call(b): [v_i] in R with i==1..nodeCount() -> scalar values # call(b): [v_i] in R³ with i==1..nodeCount() -> value values if len(test) == mesh.boundary(0).nodeCount(): if len(test[0]) == mesh.dim(): return True if test.ndim == 2 and len(test[0]) == mesh.dim(): return True except BaseException as e: if verbose: print(e) return False for key, _bVal in bc.items(): if key == 'Robin': continue if verbose: print("test vector values for:", _bVal) if isinstance(_bVal, list): if isinstance(_bVal[0], list): # [[nodeID, [x, y, z]], [nodeID, [x, y, z]]] if testForV3(_bVal[0][1]): return True else: # [nodeID, [x, y, z]] if testForV3(_bVal[1]): return True elif isinstance(_bVal, dict): # {key, bc} for key, test in _bVal.items(): if testForV3(test): return True else: # [x, y, z] = call(boundary) if testForV3(_bVal): return True return False
[docs] def parseMapToCellArray(attributeMap, mesh, default=0.0): """ Parse a value map to cell attributes. A map should consist of pairs of marker and value. A marker is an integer and corresponds to the cell.marker(). Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` For each cell of mesh a value will be returned. attributeMap : list | dict List of pairs [marker, value] ] || [[marker, value]], or dictionary with marker keys default : float [0.0] Fill all unmapped attributes to this default. Returns ------- att : array Array of length mesh.cellCount() """ # pg.warn('check if obsolete: parseMapToCellArray') att = pg.Vector(mesh.cellCount(), default) if isinstance(attributeMap, dict): for marker, value in attributeMap.items(): idx = pg.find(mesh.cellMarkers() == marker) if len(idx) == 0: pg.warn("parseMapToCellArray: cannot find marker " + str(marker) + " within mesh.") else: if isinstance(value, complex): if not isinstance(att, pg.CVector): att = pg.math.toComplex(att) att.setVal(val=value, ids=idx) else: att.setVal(val=float(value), ids=idx) elif hasattr(attributeMap, '__len__'): if not hasattr(attributeMap[0], '__len__'): # assuming [marker, value] attributeMap = [attributeMap] for pair in attributeMap: if hasattr(pair, '__len__'): idx = pg.find(mesh.cellMarkers() == pair[0]) if len(idx) == 0: pg.warn("parseMapToCellArray: cannot find marker " + str(pair[0]) + " within mesh.") else: # print('---------------------') # print(att, idx, pair[1], type(pair[1]), float(pair[1])) if isinstance(pair[1], complex): if not isinstance(att, pg.CVector): att = pg.math.toComplex(att) att.setVal(val=pair[1], ids=idx) else: att.setVal(val=float(pair[1]), ids=idx) else: raise Exception("Please provide a list of [int, value] pairs" + str(pair)) else: print("attributeMap: ", attributeMap) raise Exception("Cannot interpret attributeMap!") return att
[docs] def grad(mesh, u, r=None): r""" Return the discrete interpolated gradient :math:`\mathbf{v}` for a given scalar field :math:`\mathbf{u}`. .. math:: \mathbf{v}(\mathbf{r}_{\mathcal{C}}) &= \nabla u(\mathbf{r}_{\mathcal{N}}) \\ (\mathbf{v_x}(\mathbf{r}_{\mathcal{C}}), \mathbf{v_y}(\mathbf{r}_{\mathcal{C}}), \mathbf{v_z}(\mathbf{r}_{\mathcal{C}}))^{\text{T}} &= \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z}\right)^{\text{T}} With :math:`\mathcal{N}=\cup_{i=0}^{N} \text{Node}_i`, :math:`\mathcal{C}=\cup_{j=0}^{M} \text{Cell}_j`, :math:`\mathbf{u}=\{u(\mathbf{r}_{i})\}\in I\!R` and :math:`\mathbf{r}_{i} = (x_i, y_i, z_i)^{\text{T}}` The discrete scalar field :math:`\mathbf{u}(\mathbf{r}_{\mathcal{N}})` need to be defined for each node position :math:`\mathbf{r}_{\mathcal{N}}`. The resulting vector field :math:`\mathbf{v}(\mathbf{r}_{\mathcal{C}})` is defined for each cell center position :math:`\mathbf{r}_{\mathcal{C}}`. If you need other positions than the cell center, provide an appropriate array of coordinates :math:`\mathbf{r}`. Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` Discretization base, interpolation will be performed via finite element base shape functions. u : array | callable Scalar field per mesh node position or an appropriate callable([[x,y,z]]) r : ndarray((M, 3)) [mesh.cellCenter()] Alternative target coordinates :math:`\mathbf{r} for the resulting gradient field. i.e., the positions where the vector field is defined. Default are all cell centers. Returns ------- v : ndarray((M, 3)) Resulting vector field defined on :math:`\mathbf{v}(\mathbf{r}_{\mathcal{C}})`. M is number of cells or length of given alternative coordinates r. Examples -------- >>> import numpy as np >>> import matplotlib.pyplot as plt >>> import pygimli as pg >>> fig, ax = plt.subplots() >>> mesh = pg.createGrid(x=np.linspace(0, 1, 20), y=np.linspace(0, 1, 20)) >>> u = lambda p: pg.x(p)**2 * pg.y(p) >>> _ = pg.show(mesh, u(mesh.positions()), ax=ax) >>> _ = pg.show(mesh, [2.*pg.y(mesh.cellCenters())*pg.x(mesh.cellCenters()), ... pg.x(mesh.cellCenters())**2], ax=ax) >>> _ = pg.show(mesh, pg.solver.grad(mesh, u), ax=ax, color='w', ... linewidth=0.4) >>> plt.show() """ if r is None: r = mesh.cellCenters() uv = u if callable(u) and not isinstance(u, pg.Vector): uv = u(mesh.positions()) if len(uv) == mesh.cellCount(): uv = pg.meshtools.cellDataToNodeData(mesh, uv) v = np.ndarray((len(r), 3)) for i, _ in enumerate(v): c = mesh.findCell(r[i]) if c: v[i] = c.grad(r[i], uv) return v
[docs] def div(mesh, v): r"""Return the discrete interpolated divergence field. Return the discrete interpolated divergence field. :math:`\mathbf{u}` for each cell for a given vector field :math:`\mathbf{v}`. First order integration via boundary center. .. math:: d(cells) & = \nabla\cdot\vec{v} \\ d(c_i) & = \sum_{j=0}^{N_B}\vec{v}_{B_j} \cdot \vec{n}_{B_j} Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` Discretization base, interpolation will be performed via finite element base shape functions. V : array(N,3) | R3Vector Vector field at cell centers or boundary centers Returns ------- d : array(M) Array of divergence values for each cell in the given mesh. Examples -------- >>> import pygimli as pg >>> import numpy as np >>> v = lambda p: p >>> mesh = pg.createGrid(x=np.linspace(0, 1, 4)) >>> print(pg.math.round(pg.solver.div(mesh, v(mesh.boundaryCenters())), 1e-5)) 3 [1.0, 1.0, 1.0] >>> print(pg.math.round(pg.solver.div(mesh, v(mesh.cellCenters())), 1e-5)) 3 [0.5, 1.0, 0.5] >>> mesh = pg.createGrid(x=np.linspace(0, 1, 4), ... y=np.linspace(0, 1, 4)) >>> print(pg.math.round(pg.solver.div(mesh, v(mesh.boundaryCenters())), 1e-5)) 9 [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0] >>> divCells = pg.solver.div(mesh, v(mesh.cellCenters())) >>> # divergence from boundary values are exact where the divergence from >>> # interpolated cell center values wrong due to interpolation to boundary >>> print(sum(divCells)) 12.0 >>> mesh = pg.createGrid(x=np.linspace(0, 1, 4), ... y=np.linspace(0, 1, 4), ... z=np.linspace(0, 1, 4)) >>> print(sum(pg.solver.div(mesh, v(mesh.boundaryCenters())))) 81.0 >>> divCells = pg.solver.div(mesh, v(mesh.cellCenters())) >>> print(sum(divCells)) 54.0 """ mesh.createNeighborInfos() d = None if hasattr(v, '__len__'): if len(v) == mesh.boundaryCount(): d = mesh.divergence(v) elif len(v) == mesh.nodeCount(): d = mesh.divergence(pg.meshtools.nodeDataToBoundaryData(mesh, v)) elif len(v) == mesh.cellCount(): CtB = mesh.cellToBoundaryInterpolation() d = mesh.divergence(np.array([CtB*pg.x(v), CtB*pg.y(v), CtB*pg.z(v)]).T) else: print(len(v), mesh) raise BaseException("implement me") elif callable(v): raise BaseException("implement me") return d
[docs] def divergence(mesh, func=None, normMap=None, order=1): """Divergence for callable function func((x,y,z)). MOVE THIS to a better place Divergence for callable function func((x,y,z)). Return sum div over boundary. Parameters ---------- Returns ------- """ if func is None: func = lambda r: r diverg = 0 directionCheck = False if mesh.cellCount() > 0: directionCheck = True bNorms = None if normMap is not None: bNorms = np.zeros((mesh.boundaryCount(), 2)) for pair in normMap: bounds = mesh.findBoundaryByMarker(pair[0]) for b in bounds: bN = [0.0, 0.0] if not isinstance(pair[1][0], str): bN[0] = pair[1][0] if not isinstance(pair[1][1], str): bN[1] = pair[1][1] bNorms[b.id()] = bN for b in mesh.boundaries(): if directionCheck: if b.leftCell() is None and b.rightCell() is None: # print(b.id(), b.leftCell(), b.rightCell()) sw = pg.Stopwatch(True) mesh.createNeighborInfos() print("NeighborInfos()", sw.duration(True)) # return gauss(grid, F) # don't calc for inner boundaries if b.leftCell() is not None and b.rightCell() is not None: continue divS = 0 shape = b.shape() if order == 1: if bNorms is not None: divS = shape.norm().dot(bNorms[b.id()]) * shape.domainSize() else: divS = shape.norm().dot( func(shape.center())) * shape.domainSize() else: weights = pg.core.IntegrationRules.instance().weights(shape, order) abscissa = pg.core.IntegrationRules.instance().abscissa(shape, order) for i, p in enumerate(abscissa): rPos = shape.xyz(p) divS += shape.norm().dot(func(rPos)) * \ weights[i] * shape.domainSize() if directionCheck and b.leftCell() is None: divS *= -1 # raise Exception("invalid mesh: left is None .. every # boundary need leftCell") diverg += divS return diverg
[docs] def triDiagToeplitz(dom, a, l, r, start=0, end=-1): """Create tri-diagonal Toeplitz matrix.""" A = pg.matrix.SparseMapMatrix(dom, dom) if end == -1: end = dom for i in range(start, end): A.addVal(i, i, a) if i > start: A.addVal(i, i - 1, l) if i < end - 1: A.addVal(i, i + 1, r) return A
[docs] def identity(dom, start=0, end=-1, scale=1): """Create identity matrix.""" A = pg.matrix.SparseMapMatrix(dom, dom) if end == -1: end = dom for i in range(start, end): if hasattr(scale, '__len__'): A.addVal(i, i, scale[i]) else: A.addVal(i, i, scale) return A
[docs] def showSparseMatrix(mat, full=False): """Show the content of a sparse matrix. Parameters ---------- mat: :gimliapi:`GIMLI::SparseMatrix` | :gimliapi:`GIMLI::SparseMapMatrix` Matrix to be shown. full: bool [False] Show as dense matrix. """ if isinstance(mat, pg.matrix.SparseMapMatrix): m_ = pg.matrix.SparseMatrix(mat) return showSparseMatrix(m_, full) else: rows = mat.vecRowIdx() cols = mat.vecColPtr() vals = mat.vecVals() matD = None if full: matD = pg.Matrix(mat.rows(), mat.cols()) for i in range(mat.rows()): for j in range(cols[i], cols[i + 1]): if full: matD[i, rows[j]] = vals[j] else: if vals[j] != 0: print(i, rows[j], vals[j]) if full: print(np.array(matD))
[docs] class LinSolver(object): """Proxy class for the solution of linear systems of equations."""
[docs] def __init__(self, mat=None, solver=None, verbose=False, **kwargs): """Init the solver proxy class with a matrix and start factorization. Args ---- solver: str [None] Name for the used solver (pg (umfpack or cholmod), scipy). If solver is none decide from matrix type. """ self._m = None # hold local copy if we need to convert the matrix self.verbose = verbose self._solver = None self.factorTime = 0.0 self.solvingTime = 0.0 self.solver = '' self._factorize = 'factorizePG' self._factorized = False self._desiredArrayType = np.array if solver is None: if isinstance(mat, pg.matrix.MatrixBase): solver = 'PG' elif isinstance(mat, np.ndarray): solver = 'numpy' pg.critical("Not yet implemented!") else: from scipy.sparse import spmatrix if isinstance(mat, spmatrix): solver = 'SciPy' if solver.lower() == 'pg': self.solver = 'PG' elif solver.lower() == 'scipy': self.solver = 'SciPy' else: self.solver = solver self._factorize = 'factorize' + self.solver if self.verbose: pg.info("Solving with {0}".format(self.solver)) if mat is not None: self.factorize(mat)
[docs] def isFactorized(self): return self._factorized
[docs] def factorize(self, mat): swatch = pg.Stopwatch() getattr(self, self._factorize)(mat) self.factorTime = swatch.duration(restart=True) if self.verbose: pg.info("Matrix factorization:", self.factorTime) self._factorized = True
[docs] def factorizePG(self, mat): """""" self._m = pg.utils.toSparseMatrix(mat) self._desiredArrayType = pg.Vector self._solver = pg.core.LinSolver(self._m, verbose=self.verbose)
[docs] def factorizeSciPy(self, mat): """""" self._m = pg.utils.sparseMatrix2csr(mat) # scipy is not dependency # scipy = pg.optImport('scipy', 'Used for sparse linear solver.') from scipy.sparse.linalg import factorized self._desiredArrayType = np.array self._solver = factorized(self._m)
def __call__(self, b): """short cut to self.solve(b)""" return self.solve(b) def _convertRHS(self, b): """Convert right hand side vector into the desired format.""" if not isinstance(b, type(self._desiredArrayType(0))): return self._desiredArrayType(b) return b
[docs] def solve(self, b): """ """ swatch = pg.Stopwatch() x = self._solver(self._convertRHS(b)) self.solverTime = swatch.duration(restart=True) if self.verbose: pg.info("Matrix solve:", self.solverTime) return x
[docs] def linSolve(mat, b, solver=None, verbose=False, **kwargs): r"""Direct linear solution after :math:`\textbf{x}` using core LinSolver. .. math:: \textbf{A}\textbf{x} = \textbf{b} If :math:`\textbf{A}` is symmetric, sparse and positive definite. Parameters ---------- mat: :gimliapi:`GIMLI::RSparseMatrix`, :gimliapi:`GIMLI::RSparseMapMatrix`, numpy.array System matrix. Need to be symmetric, sparse and positive definite. b: iterable array Right hand side of the equation. solver: str [None] Try to choose a solver, 'pg' for pygimli core cholmod or umfpack. 'np' for numpy linalg or scipy.sparse.linalg. Automatic choosing if solver is None depending on matrixtype. verbose: bool [False] Be verbose. Returns ------- x: :gimliapi:`GIMLI::Vector` Solution vector. """ # TODO!! refactor with LinSolver swatch = pg.Stopwatch() reorder = kwargs.pop('reorder', False) # perm = None # determine the solver if none set if solver is None: if isinstance(mat, pg.matrix.MatrixBase): solver = 'pg' elif isinstance(mat, np.ndarray): solver = 'numpy' else: from scipy.sparse import spmatrix if isinstance(mat, spmatrix): solver = 'scipy' if solver == 'pg': # core proxy to cholmod and LDL for float and umfpack for complex if reorder is True: pg.warning( 'Matrix reordering for pg core solver not yet implemented') _m = pg.utils.toSparseMatrix(mat) solver = pg.core.LinSolver(_m, verbose=verbose) if verbose: pg.info("Solving with {0}".format(solver.solverName())) pg.info("Matrix factorization:", swatch.duration(restart=True)) x = solver.solve(b) if verbose: pg.info("Matrix solution:", swatch.duration()) elif solver == 'numpy': if verbose: pg.info("Solving with np.linalg.solve") x = np.linalg.solve(mat, b) elif solver == 'scipy': # pg._r(swatch.duration(restart=True)) _m = pg.utils.sparseMatrix2csr(mat) # pg._r('convert', swatch.duration(restart=True)) # scipy is now a dependency # scipy = pg.optImport('scipy', 'Used for sparse linear solver.') # pg._r('import', swatch.duration(restart=True)) from scipy.sparse.linalg import spsolve if verbose: pg.info("Solving with scipy.sparse.spsolve") if reorder is True and 0: def permCOO(M, perm): # M.indices = perm.take(M.indices) # M = M.tocsc() # M.indices = perm.take(M.indices) # return M.tocsr() return M[np.ix_(perm, perm)] # print(M.row.shape, M.col.shape) # rowP = perm[M.row] # colP = perm[M.col] # print(rowP.shape, colP.shape) # MP = scipy.sparse.coo_matrix((M.data, (rowP, colP)), # shape=M.shape) # return MP # perm = scipy.sparse.csgraph.reverse_cuthill_mckee(_m) # pg._r('reverse_cuthill_mckee', swatch.duration(restart=True)) # ax, _ = pg.show(_m) # _m = permCOO(_m, perm) # _m = permCOO(_m.tocoo(), perm).tocsr() # pg.show(_m, ax=ax, color='green') # pg._r('perm matrix', swatch.duration(restart=True)) # x = spsolve(_m, b.array())[perm] # x = spsolve(_m, b.array()[perm])#[perm] # x = x[perm] else: x = spsolve(_m, b) # pg._r(swatch.duration()) return x
[docs] def applyDirichlet(mat, rhs, uDirIndex, uDirichlet): """This should be moved directly into the core""" # idx = np.argsort(uDirIndex) # print(idx) # print(np.asarray(uDirIndex)[idx]) # print(np.asarray(uDirichlet[idx])) if mat is not None: if rhs is not None: uDir = pg.Vector(mat.rows(), 0.0) uDir.setVal(uDirichlet, uDirIndex) rhs -= mat * uDir for i in uDirIndex: mat.cleanRow(i) mat.cleanCol(i) mat.setVal(i, i, 1.0) if rhs is not None: rhs[uDirIndex] = uDirichlet
# rhs.setVal(uDirichlet, uDirIndex)
[docs] def getDirichletMap(mat, boundaryPairs, time=0.0, userData={}, nodePairs=None, dofOffset=0, nCoeff=1, dofPerCoeff=None): r"""Get map of index: dirichlet value Apply Dirichlet boundary condition to the system matrix S and rhs vector. The right hand side values for h can be given for each boundary element individually by setting proper boundary pair arguments. .. math:: u(\textbf{r}, t) = h \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega= \Gamma_{\text{Dirichlet}} Parameters ---------- mat: :gimliapi:`GIMLI::RSparseMatrix` System matrix of the system equation. boundaryPair: list() List of pairs [:gimliapi:`GIMLI::Boundary`, h]. The value :math:`h` will assigned to the nodes of the boundaries. Later assignment overwrites prior. :math:`h` need to be a scalar value (float or int) or a value generator function that will be executed at runtime. See :py:mod:`pygimli.solver.solver.parseArgToBoundaries` and :ref:`tut:modelling_bc` for example syntax, nodePairs: list() | callable List of pairs [nodeID, uD]. The value uD will assigned to the nodes given there ids. This node value settings will overwrite any prior settings due to boundaryPair. time: float Will be forwarded to value generator. userData: class Will be forwarded to value generator. dofOffset: int[0] Offset for matrix index. """ if not hasattr(boundaryPairs, '__getitem__'): raise BaseException("Boundary pairs need to be a list of " "[boundary, value]") # uDirNodes = [] # [] uDirVal = dict() # {nID: val} def _genVecUd(n, ud, dofOffset, nCoeff=1, dofPerCoeff=None): ret = {} if callable(ud): pg.error("callable node pairs need to be implemented.") if isinstance(n, pg.core.Node): idx = dofOffset + n.id() else: idx = dofOffset + n if hasattr(ud, '__iter__'): # vector valued problem if dofPerCoeff is None: if mat.shape[0] % len(ud) != 0: print(mat) print(mat.shape, len(ud)) pg.error("Matrix size missmatch for vector valued problem") else: dofPerCoeff = mat.shape[0] // len(ud) if nCoeff == 1: nCoeff = len(ud) for i in range(nCoeff): if ud[i] is not None: ret[idx + i * dofPerCoeff] = ud[i] else: if nCoeff > 1: print('nCoeff:', nCoeff, 'ud:', ud, 'idx:', idx) pg.error('number of coefficents > 1 but uDirichlet is scalar.') if ud is not None: ret[idx] = ud return ret for pair in boundaryPairs: ent = pair[0] val = pair[1] # print('**', ent, val) uD = generateBoundaryValue(ent, val, time, userData, nCoeff=nCoeff) # print('\t', uD) if uD is not None: if isinstance(ent, pg.core.Node): uDirVal.update(_genVecUd(ent, uD, dofOffset)) else: if isinstance(uD, float): pg.critical(uD) uD = [uD] * ent.nodeCount() if len(uD) == ent.nodeCount(): # print('uD', uD, nCoeff, dofPerCoeff) for i, n in enumerate(ent.nodes()): uDirVal.update(_genVecUd(n, uD[i], dofOffset, nCoeff=nCoeff, dofPerCoeff=dofPerCoeff)) else: pg.error('Dirichlet values per boundary need to have ' 'length of boundary.nodeCount()') if nodePairs is not None: # print("nodePairs", nodePairs) if len(nodePairs) == 2 and isinstance(nodePairs[0], int): # assume a single Node [NodeId, val] nodePairs = [nodePairs] for [n, val] in nodePairs: uDirVal.update(_genVecUd(n, val, dofOffset, nCoeff=nCoeff, dofPerCoeff=dofPerCoeff)) return uDirVal
[docs] def assembleDirichletBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, nodePairs=None, dofOffset=0, nCoeff=1, dofPerCoeff=None): r"""Apply Dirichlet boundary condition. Args ---- rhs: :py:mod:`Vector` Right hand side vector of the system equation will bet set to :math:`u_{\text{D}}` """ uDirVal = getDirichletMap(mat, boundaryPairs, time=time, userData=userData, nodePairs=nodePairs, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=dofPerCoeff) # pg._g(list(uDirVal.keys()), list(uDirVal.values())) if not uDirVal.keys(): return applyDirichlet(mat, rhs, list(uDirVal.keys()), list(uDirVal.values())) return uDirVal
[docs] def assembleNeumannBC(rhs, boundaryPairs, nDim=1, time=0.0, userData={}, dofOffset=0, nCoeff=1, dofPerCoeff=None): r"""Apply Neumann condition to the system matrix S. Apply Neumann condition to the system matrix S. The right hand side values for g can be given for each boundary element individually by setting proper boundary pair arguments. .. math:: \frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} = \textbf{n}\nabla u(\textbf{r}, t) = g \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Neumann}} Parameters ---------- rhs: :py:mod:`Vector` Right hand side vector of length node count. boundaryPair : list() List of pairs [ :gimliapi:`GIMLI::Boundary`, g ]. The value :math:`g` will assigned to the nodes of the boundaries. Later assignment overwrites prior. :math:`g` need to be a scalar value (float or int) or a value generator function that will be executed at run time. See :py:mod:`pygimli.solver.solver.parseArgToBoundaries` and :ref:`tut:modelling_bc` for example syntax, nDim: int [1] Number of dimensions for vector valued problems. The rhs array need to have the correct size, i.e., number of Nodes * mesh.dimension() time: float Will be forwarded to value generator. userData: class Will be forwarded to value generator. dofOffset: int[0] Offset for matrix index. """ if rhs is None: raise BaseException("Neumann Boundary condition needs rhs vector.") if not hasattr(boundaryPairs, '__getitem__'): raise BaseException("Boundary pairs need to be a list of " "[boundary, value]") Se = pg.matrix.ElementMatrix() dof = len(rhs) // nDim for pair in boundaryPairs: boundary = pair[0] val = pair[1] # print('+++++', boundary) # print('\t', val) g = generateBoundaryValue(boundary, val, time, userData, nCoeff=nCoeff) # print('\t', g) # if a is not None: # pg.warning('Scaling of Neumann values necessary? Check!') # try: # g *= a[boundary.leftCell().id()] # except BaseException as e: # print(boundary.leftCell()) # print(boundary.leftCell().id()) # print(len(a)) # pg.warn('Insufficient cell information.') if g is not None: Se.u(boundary) for dim in range(nDim): if nDim == 1: gd = g else: if isinstance(g, list) and len(g) == nDim: gd = g[dim] else: gd = g.T[dim] idx = Se.ids() + dim*dof + dofOffset if isinstance(gd, float) and gd == 0: continue if hasattr(gd, '__iter__') and not np.any(gd): continue # print(nDim, g, gd) if isinstance(rhs, pg.Vector): # print(Se) # pg.info(sum(Se.row(0))) # pg.info(Se.row(0), gd, idx) rhs.addVal(Se.row(0) * gd, idx) # rhs.setVal(Se.row(0) * gd, idx) # rhs.add(Se, g) else: # check # pg.error('check') rhs[idx] += Se.row(0) * gd
# for i, j in enumerate(Se.ids()): # rhs[j + dim*dof] += Se.row(0)[i] * gd
[docs] def assembleRobinBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, dofOffset=0, nCoeff=1, dofPerCoeff=None): r"""Apply Robin boundary condition. Apply Robin boundary condition to the system matrix and the rhs vector: .. math:: \frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} & = \alpha(u_0-u) \quad\text{or} \\ \beta\frac{\partial u(\textbf{r}, t)}{\partial\textbf{n}} + \alpha u & = \gamma \\ & \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega= \Gamma_{\text{Robin}}\\ Parameters ---------- mat: :gimliapi:`GIMLI::SparseMatrix` System matrix of the system equation. boundaryPair: list List of pairs [:gimliapi:`GIMLI::Boundary`, :math:`a, u_0` | :math:`\alpha, \beta, \gamma`]. The values will assigned to the nodes of the boundaries. Later assignment overwrites prior. Values can be a single value for :math:`\alpha` or :math:`a`, two values will be interpreted as :math:`a, u_0`, and three values will be :math:`\alpha, \beta, \gamma`. Also generator (callable) is possible which will be executed at runtime See :py:mod:`pygimli.solver.solver.parseArgToBoundaries` :ref:`tut:modelling_bc` or testing/test_FEM.py for example syntax. time: float Will be forwarded to value generator. userData: dict Will be forwarded to value generator. dofOffset: int[0] Offset for matrix index. """ if not hasattr(boundaryPairs, '__getitem__'): raise BaseException("Boundary pairs need to be a list of " "[boundary, value]") S_Dir = pg.matrix.ElementMatrix() S_Neu = pg.matrix.ElementMatrix() # if isinstance(rhs, np.ndarray): # rhs = pg.Vector(rhs) for pair in boundaryPairs: boundary = pair[0] val = pair[1] # print('val:', val) # du/dn = a(u0-u) || \beta du/dn + \alpha u = \gamma # combines to Matrix + au = RHS + au0 u0 = None a = generateBoundaryValue(boundary, val, time, userData, expectList=True, nCoeff=nCoeff) try: if a.ndim == 2 and len(a) == boundary.nodeCount(): a = a[0] except BaseException: # expecting [[a| a, u0 | a b g]_i] for i in boundary.nodes() print(boundary) print(a) print(a.ndim) pg.error("Can't interprete robin value.") if hasattr(a, '__iter__'): if len(a) == 1: a = a[0] elif len(a) == 2: u0 = a[1] a = a[0] elif len(a) == 3: alpha, beta, gamma = a[0], a[1], a[2] # a = [alpha, beta, gamma] if alpha != 0: u0 = gamma/alpha else: pg.warn('Robin boundary condition parmeter alpha is zero, ' 'falling back to Neumann condition.') u0 = 0.0 if beta != 0: a = alpha/beta else: pg.warn('Robin boundary condition parmeter beta is zero, ' 'please consider using Dirichlet instead.') a = 0.0 if a is not None and a != 0.0: S_Dir.u2(boundary) mat.add(S_Dir, scale=a) # Sp *= p # S += Sp if u0 is not None and u0 != 0.0: S_Neu.u(boundary) rhs.add(S_Neu, a * u0)
[docs] def assembleBC(bc, mesh, mat, rhs, time=None, userData={}, dofOffset=0, nCoeff=1): r"""Shortcut to apply all boundary conditions. Shortcut to apply all boundary conditions will only forward to appropriate assemble functions. Parameters ---------- Returns ------- map{id: uDirichlet}: Map of index to Dirichlet value. None """ # we can't iterate because we want the following fixed order dirichletMap = {} bct = dict(bc) nDim = 1 if mat is not None: if mat.rows() == mesh.nodeCount() * mesh.dim(): nDim = mesh.dim() if 'Neumann' in bct: assembleNeumannBC(rhs, parseArgToBoundaries(bct.pop('Neumann'), mesh), nDim=nDim, time=time, userData=userData, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=mesh.nodeCount()) if 'Robin' in bct: assembleRobinBC(mat, parseArgToBoundaries(bct.pop('Robin'), mesh), rhs=rhs, time=time, userData=userData, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=mesh.nodeCount()) if 'Dirichlet' in bct: uD = assembleDirichletBC( mat, parseArgToBoundaries(bct.pop('Dirichlet'), mesh), rhs=rhs, time=time, userData=userData, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=mesh.nodeCount()) dirichletMap.update(uD) if 'Nodes' in bct: # 'Nodes' : [list(Nodes), callable(Node)] ## for selected Nodes # 'Nodes' : callable(Node) ## for all nodes bc = bct.pop('Nodes') if isinstance(bc, list): nodes = bc[0] val = bc[1] else: nodes = mesh.nodes() val = bc nP = [] if callable(val): for n in nodes: nP.append([n.id(), val(n)]) else: pg.critical("Nodes boundary need a callable(Node)") uD = assembleDirichletBC( mat, [], nodePairs=nP, rhs=rhs, time=time, userData=userData, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=mesh.nodeCount()) dirichletMap.update(uD) if 'Node' in bct: uD = assembleDirichletBC( mat, [], nodePairs=bct.pop('Node'), rhs=rhs, time=time, userData=userData, dofOffset=dofOffset, nCoeff=nCoeff, dofPerCoeff=mesh.nodeCount()) dirichletMap.update(uD) if len(bct.keys()) > 0: pg.warn("Unknown boundary condition[s]" + str(bct.keys()) + " will be ignored") return dirichletMap
[docs] def assembleLoadVector(mesh, f, userData={}): r"""Assemble the load vector. See createLoadVector.""" pg.deprecate('createLoadVector') # 20200115 return createLoadVector(mesh, f, userData)
[docs] def createForceVector(mesh, f, userData={}): """ Create a right hand side vector for vector valued solutions. Parameters ---------- f: [ convertable ] List of rhs side options. Must be convertable to createLoadVector. See :py:mod:`createLoadVector` rhs: np.array() Squeezed vector of length mesh.nodeCount() * mesh.dimensions() """ if not isinstance(f, list): pg.error("Create Force Vector need list of attribute f with an entry " "for each dimension.") rhs = np.zeros(mesh.nodeCount() * mesh.dim()) for i in range(mesh.dim()): rhs[i*mesh.nodeCount():(i+1)*mesh.nodeCount()] = \ createLoadVector(mesh, f[i], userData) # rhs.reshape(mesh.nodeCount() * mesh.dim()) #contiguity not guarantied return rhs
[docs] def createLoadVector(mesh, f=1.0, userData={}): """Create right hand side vector based on the given mesh and load values (scalar solution) or force vectors (vector value solution). Create right hand side based on the given mesh and load or force values. TODO ---- * Callable for vector valued problems * Callable called dynamic on demand Parameters ---------- f: float[1.0], array, callable(cell, [userData]), [f_x, f_y, f_z] * float will be assumed as constant for all cells like rhs = rhs(np.ones(mesh.cellCount() * f), * array of length mesh.cellCount() will be processed as load value for each cell: rhs = rhs(f), * array of length mesh.nodeCount() is assumed to be already processed correct: rhs = f * callable is evaluated on once for each cell and need to return a load value for each cell and can have optional a userData dictionary: `f_cell = f(cell, [userData={}])` rhs = rhs(f(c, userData) for c in mesh.cells()) * list with length of mesh.dimension() of float or array entries will create a squeezed rhs for vector valued problems rhs = squeeze([rhs(f[0]), rhs(f[1]), rhs(f[2])]) Returns ------- rhs: pg.Vector(mesh.nodeCount()) Right-hand side load vector for scalar values or squeezed vector values """ # f is dict('Node':callable, 'Cell': callable) if isinstance(f, dict): if 'Node' in f: fn = [] if callable(f['Node']): for n in mesh.nodes(): fn.append(f['Node'](n, **userData)) if hasattr(fn[0], '__iter__'): # result is vector valued return createLoadVector(mesh, [fi for fi in np.array(fn).T], userData=userData) return createLoadVector(mesh, fn, userData=userData) elif 'Cell' in f: pg.error('Implement me!, createLoadVector()') # fix for the lazy if isinstance(f, int): f = float(f) # f is list [fx, fy, [fz]] for vector problems if isinstance(f, list): if len(f) == mesh.dim(): return createForceVector(mesh, f, userData=userData) # f is list of array [f_0, f_1, ..., f_n] for scalar problems if isinstance(f, list) or hasattr(f, 'ndim'): if isinstance(f, list): rhs = np.zeros((len(f), mesh.nodeCount())) for i, fi in enumerate(f): userData['i'] = i rhs[i] = createLoadVector(mesh, fi, userData) return rhs elif f.ndim == 2: # assume rhs [n, nNodes] array is already a valid if len(f[0]) == mesh.nodeCount(): return f rhs = pg.Vector(mesh.nodeCount(), 0) fArray = None if hasattr(f, '__len__'): if len(f) == mesh.cellCount(): # scalar values for each cell fArray = f elif len(f) == mesh.nodeCount(): # scalar values for each node fArray = f elif len(f) == mesh.nodeCount() * mesh.dim(): # vector values for each node # maybe just for special cases with allready processed rhs return f elif callable(f) and not isinstance(f, pg.Vector): fArray = pg.Vector(mesh.cellCount()) for c in mesh.cells(): if userData is not None and userData.keys(): fArray[c.id()] = f(c, userData) else: fArray[c.id()] = f(c) if fArray is None: fArray = cellValues(mesh, f, userData=userData) if len(fArray) == mesh.cellCount(): b_l = pg.matrix.ElementMatrix() for c in mesh.cells(): if fArray[c.id()] != 0.0: b_l.u(c) rhs.add(b_l, fArray[c.id()]) # print("test reference solution:") # rhsRef = pg.Vector(mesh.nodeCount(), 0) # for c in mesh.cells(): # b_l.u(c) # for i, idx in enumerate(b_l.idx()): # rhsRef[idx] += b_l.row(0)[i] * fArray[c.id()] # np.testing.assert_allclose(rhs, rhsRef) # print("Remove revtest in assembleLoadVector after check") elif len(fArray) == mesh.nodeCount(): # nodal values fA = pg.Vector(fArray) b_l = pg.matrix.ElementMatrix() for c in mesh.cells(): b_l.u(c) # rhs.addVal(b_l.row(0) * fArray[b_l.idx()], b_l.idx()) rhs.add(b_l, fA) # print("test reference solution:") # rhsRef = pg.Vector(mesh.nodeCount(), 0) # for c in mesh.cells(): # b_l.u(c) # for i, idx in enumerate(b_l.idx()): # rhsRef[idx] += b_l.row(0)[i] * fA[idx] # np.testing.assert_allclose(rhs, rhsRef) # print("Remove revtest in assembleLoadVector after check", # sum(rhs), sum(rhsRef)) # rhs = pg.Vector(fArray) else: raise Exception("Load vector have the wrong size: " + str(len(fArray))) return rhs
[docs] def createStiffnessMatrix(mesh, a=None, isVector=False): r"""Create the Stiffness matrix. Calculates the Stiffness matrix :math:`{\bf S}` for the given mesh scaled with the per cell values a. ..math:: ... Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` Arbitrary mesh to calculate the stiffness for. Type of base and shape functions depends on the cell types. a : iterable of type float, int, complex, RMatrix, CMatrix Per cell values., e.g., physical parameter. Length of a need to be mesh.cellCount(). If None given default is 1. isVector : bool [False] We want to solve for vector valued problems. Resulting SparseMatrix is a SparseMapMatrix and have the dimension (nNodes * nDims, nNodes * nDims) with nNodes = mesh.nodeCount() and nDims = mesh.dimension(). Returns ------- A : :gimliapi:`GIMLI::[C]SparseMatrix` | [C]SparseMapMatrix Stiffness matrix, with real or complex values. """ if mesh.cellCount() == 0: print(mesh) raise Exception("Mesh invalid") if a is None: a = pg.Vector(mesh.cellCount(), 1.0) A = None if isVector is False: if isinstance(a[0], float) or \ isinstance(a[0], int) or \ isinstance(a[0], np.float64): A = pg.matrix.SparseMatrix() A.fillStiffnessMatrix(mesh, a) return A dof = 0 nDof = mesh.nodeCount() else: dof = mesh.nodeCount() nDof = mesh.nodeCount() * mesh.dimension() # if vector or scalar(Complex) if pg.isComplex(a[0]): isComplex = True A = pg.matrix.CSparseMapMatrix(nDof, nDof) else: isComplex = False A = pg.matrix.SparseMapMatrix(nDof, nDof) al = pg.core.ElementMatrix(dof=dof) if len(a) != mesh.cellCount(): pg.error('Number of cell values need to match cell count') for c in mesh.cells(): if isComplex is True: # al.gradU2(c, 1.0) al.ux2uy2uz2(c) A.add(al, scale=a[c.id()]) else: if pg.isScalar(a[c.id()]): al.gradU2(c, a[c.id()]) A.add(al) else: if hasattr(a[c.id()], 'voigtNotation'): vN = a[c.id()].voigtNotation else: vN = False # al.gradU2(c, a[c.id()], voigtNotation=vN) al.gradU2(c, np.array(a[c.id()]), voigtNotation=vN) A.add(al) if isComplex is True: return pg.matrix.CSparseMatrix(A) return pg.matrix.SparseMatrix(A)
[docs] def createMassMatrix(mesh, b=None): r"""Create the mass matrix. Calculates the Mass matrix (Finite element identity matrix) the given mesh. ..math:: ... Parameters ---------- mesh : :gimliapi:`GIMLI::Mesh` Arbitrary mesh to calculate the mass element matrix. Type of base and shape functions depends on the cell types. b : array Per cell values. If None given default is 1. Returns ------- A : :gimliapi:`GIMLI::RSparseMatrix` Mass element matrix """ # need callable here if b is None: b = pg.Vector(mesh.cellCount(), 1.0) elif not hasattr(b, '__iter__'): b = pg.Vector(mesh.cellCount(), b) B = pg.matrix.SparseMatrix() B.fillMassMatrix(mesh, b) return B
# create matrix structure regarding the mesh # B.buildSparsityPattern(mesh) # define a local element matrix # B_l = pg.matrix.ElementMatrix() # for c in mesh.cells(): # B_l.u2(c) # # check if b[i] == B*b # B_l *= b[c.id()] # B += B_l # return B
[docs] def intDomain(u, mesh=None): r"""Return integral over nodal solution :math:`u`. .. math:: \int_{\Omega} u TODO ---- * refactor * better name? * Documentation """ if mesh is not None: r = createLoadVector(mesh) return sum(r*u) pg.critical('Need a mesh to calculate the integral over domain')
def _feNorm(u, mat): """Create a norm within a Finite Element space. Create the Finite Element Norm with a preassembled system matrix. """ return np.sqrt(pg.math.dot(u, mat.mult(u)))
[docs] def normL2(u, mat=None, mesh=None): r"""Create Lebesgue (L2) norm for finite element space. Find the L2 Norm for a solution for the finite element space. :math:`u` exact solution :math:`{\bf M}` Mass matrix, i.e., Finite element identity matrix. .. math:: L2(f(x)) = || f(x) ||_{L^2} & = (\int |f(x)|^2 \mathrm{d}\:x)^{1/2} \\ & \approx h (\sum |f(x)|^2 )^{1/2} \\ L2(u) = || u ||_{L^2} & = (\int |u|^2 \mathrm{d}\:x)^{1/2} \\ & \approx (\sum M (u))^{1/2} \\ e_{L2_rel} = \frac{L2(u)}{L2(u)} & = \frac{(\sum M(u))^{1/2}}{(\sum M u)^{1/2}} The error for any approximated solution :math:`u_h` correlates to the L2 norm of 'L2Norm(u - u_h, M)'. If you like relative values, you can also normalize this error with 'L2Norm(u - u_h, M) / L2Norm(u, M)*100'. Parameters ---------- u : iterable Node based value to compute the L2 norm for. mat : Matrix Mass element matrix. mesh : :gimliapi:`GIMLI::Mesh` Mesh with the FE space to generate M if necessary. Returns ------- ret : float :math:`L2(u)` norm. """ if isinstance(mat, pg.Mesh): mesh = mat mat = None if mat is None and mesh is not None: mat = createMassMatrix(mesh) if mat is None: pg.warning("No Stiffness matrix or a mesh here, to calculate L2-Norm. " "Returning algebraic l2.") # M is Identity matrix return np.sqrt(pg.math.dot(u, u)) return _feNorm(u, mat)
[docs] def normH1(u, mat=None, mesh=None): r"""Create (H1) norm for the finite element space. Parameters ---------- u : iterable Node based value to compute the H1 norm for. mat : Matrix Stiffness matrix. mesh : :gimliapi:`GIMLI::Mesh` Mesh with the FE space to generate S if necessary. Returns ------- ret : float :math:`H1(u)` norm. """ if isinstance(mat, pg.Mesh): mesh = mat mat = None if mat is None and mesh is not None: mat = pg.solver.createStiffnessMatrix(mesh) if mat is None: raise Exception("Need Stiffness matrix or mesh to calculate H1 norm") return _feNorm(u, mat)
[docs] def solve(mesh, **kwargs): r"""Solve partial differential equation. This is a syntactic sugar convenience function for solving partial differential equation on a given mesh. Using the best guess method for the given parameter. Currently only Finite Element calculation using :py:mod:`pygimli.solver.solveFiniteElements` TODO :py:mod:`pygimli.solver.solveFiniteVolume` """ return solveFiniteElements(mesh, **kwargs)
[docs] def solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None, times=None, c=1.0, userData={}, verbose=False, **kwargs): r"""Solve partial differential equation with Finite Elements. This is a syntactic sugar convenience function for using the Finite Element functionality of the library core to solve partial differential equation (PDE) that match the following form: .. math:: c \frac{\partial u}{\partial t} & = \nabla\cdot(a \nabla u) + b u + f(\mathbf{r},t)~~|~~\Omega_{\text{Mesh}}\\ u & = h~~|~~\Gamma_{\text{Dirichlet}}\\ a\frac{\partial u}{\partial \mathbf{n}} & = g~~|~~\Gamma_{\text{Neumann}}\\ \alpha u + \beta\frac{\partial u}{\partial \mathbf{n}} & = \gamma~~|~~\Gamma_{\text{Robin}}\\ \frac{\partial u}{\partial \mathbf{n}} & = \alpha(u_0-u)~~|~~\Gamma_{\text{Robin}} for the scalar :math:`u(\mathbf{r}, t)` or vector :math:`\mathbf(u)(\mathbf{r}, t)` solution at each node of a given mesh. The Domain :math:`\Omega` and the Boundary :math:`\Gamma` are defined through the mesh with appropriate boundary marker. To ensure vector solution, either set vector forces or at least one vector component boundary condition. TODO ---- * unsteady ub and dub * 'Infinity' Boundary condition (u vanishes at infinity) * 'Cauchy' Boundary condition (guaranties u and du on same boundary) will never work here because the problem becomes ill posed and would need some inverse strategy to solve. * Example for * elastic parameter * anisotropic (float/complex) * dynamic boundary conditions * dynamic load vector * nonlinearity Parameters ---------- mesh: :gimliapi:`GIMLI::Mesh` Mesh represents spatial discretization of the calculation domain a: value | array | callable(cell, userData) Cell values of type float or complex can be scalar, anisotropy matrix or elastic tensor. b: value | array | callable(cell, userData) [None] Cell values. None means the term is unused. c: value | array | callable(cell, userData) [None] Scale the unsteady term, only for times is not None. f: value | array(cells) | array(nodes) | callable(args, kwargs) force values, for vector fields use (n x dim) values. bc: dict() Dictionary of boundary conditions. Current supported boundary conditions by dictionary keys: 'Dirichlet', 'Neumann', 'Robin', 'Node'. The dictionary can contain multiple "key: Arg" Arg will be parsed by :py:mod:`pygimli.solver.solver.parseArgPairToBoundaryArray` If the dictionary key is 'Node' then fixed values for single node indices can by be given. e.g., bc={'Node': [nodeID, value]}. Note this is only a shortcut for bc={'Dirichlet': [mesh.node(nodeID), value]}. The parameter $a$ for Neumann boundary condition is choosen automatically from the diffusivity parameter $a$ of the associated cell. times: array [None] Solve as time dependent problem for the given times. Keyword Arguments ----------------- **kwargs u0: value | array | callable(pos, userData) Node values theta: float [1] * :math:`theta = 0` means explicit Euler, maybe stable for :math:`\Delta t \quad\text{near}\quad h` * :math:`theta = 0.5`, Crank-Nicolson scheme, maybe instable * :math:`theta = 2/3`, Galerkin scheme * :math:`theta = 1`, implicit Euler If unsure choose :math:`\theta = 0.5 + \epsilon` (probably stable). dynamic: bool [False] Boundary conditions for time depending problems will be considered dynamic for each time step. stats: bool Give some statistics. progress: bool Give some calculation progress. assembleOnly: bool Stops after matrix asssemblation. Returns the system matrix A and the rhs vector. fixPureNeumann: bool [auto] If set or detected automatic, we add the additional condition: :math:`\int_\Omega u dv = 0` making elliptic problems well-posed. rhs: iterable Pre assembled rhs. Will preferred on any f settings. ws: dict The WorkSpace is a dictionary that will get some temporary data during the calculation. Any keyvalue 'u' in the dictionary is used for the resulting array. vectorValued: bool (False) Solution forced to vector valued, in case the auto detection fails Returns ------- u: array Returns the solution u either 1,n array for stationary problems or for m,n array for m time steps See also -------- :ref:`tut:modelling` and :py:mod:`pygimli.solver.solve` Examples -------- >>> # no need to import matplotlib, pygimli show does. >>> import pygimli as pg >>> import pygimli.meshtools as mt >>> world = mt.createWorld(start=[-10, 0], end=[10, -10], ... marker=1, worldMarker=False) >>> c1 = mt.createCircle(pos=[0.0, -5.0], radius=3.0, area=.1, marker=2) >>> mesh = mt.createMesh([world, c1], quality=34.3) >>> u = pg.solver.solveFiniteElements(mesh, a={1: 100.0, 2: 1.0}, ... bc={'Dirichlet':{4: 1.0, 3: 0.0}}) >>> ax = pg.show(mesh, u, showMesh=True)[0] >>> _ = pg.show(c1, ax=ax, fillRegion=False) """ if bc is None: bc = {} workSpace = kwargs.pop('ws', dict()) debug = kwargs.pop('debug', False) stats = kwargs.pop('stats', False) mesh.createNeighborInfos() if verbose: print("Mesh: ", str(mesh)) # scalar solution default vectorValues = False dof = mesh.nodeCount() # check if force vector is a vector rhs = kwargs.pop('rhs', createLoadVector(mesh, f, userData=userData)) # pg._g('###############') if len(rhs) > dof or kwargs.pop('vectorValued', _bcIsForVectorValues(bc, mesh)): if verbose: print("Solve vector valued.") vectorValues = True dof = mesh.nodeCount() * mesh.dimension() # pg._g('###############', dof) rhs.resize(dof) swatch = pg.core.Stopwatch(True) # check for material parameter # a = parseArgToArray(a, nDof=mesh.cellCount(), mesh=mesh, userData=userData) a = cellValues(mesh, a, userData=userData) isComplex = False if pg.utils.isComplex(a): isComplex = True rhs = np.array(rhs, dtype=complex) S = createStiffnessMatrix(mesh, a, isVector=vectorValues) M = None if b is not None and b != 0: b = cellValues(mesh, b, userData=userData) M = createMassMatrix(mesh, b) # pg.warn("check me") A = S - M else: A = S if times is None: if len(list(bc.items())) == 0 or \ (len(list(bc.items())) == 1 and list(bc.keys())[0] == 'Neumann'): pn = True else: pn = False fixPureNeumann = kwargs.pop('fixPureNeumann', pn) assembleBC(bc, mesh, A, rhs, time=None, userData=userData) u = None if 'u' in workSpace: u = workSpace['u'] singleForce = True if hasattr(rhs, 'ndim'): if rhs.ndim == 2: singleForce = False if u is None: u = np.zeros(rhs.shape) else: if isinstance(a[0], complex): if u is None: u = pg.CVector(rhs.size(), 0.0) rhs = pg.math.toComplex(rhs) else: if u is None: u = pg.Vector(rhs.size(), 0.0) if fixPureNeumann is True: pg.info('Fixing pure Neumann boundary condition by forcing: ' 'intDomain(u, mesh) = 0') r = createLoadVector(mesh) A = pg.BlockMatrix() A.add(S, 0, 0) A.add(r, 0, mesh.nodeCount()) A.add(r, mesh.nodeCount(), 0, transpose=True) rhs = pg.cat(rhs, pg.Vector(1, 0)) assembleTime = swatch.duration(True) if stats: stats.assembleTime = assembleTime if verbose: print("Assembling time: ", assembleTime) workSpace['Stiffness matrix'] = S workSpace['Mass matrix'] = M workSpace['System matrix'] = A workSpace['rhs'] = rhs if 'assembleOnly' in kwargs: return A, rhs if fixPureNeumann: if singleForce: uc = pg.solver.linSolve(A, rhs, 'scipy') u = uc[0:mesh.nodeCount()] else: pg.critical( 'Non-single force for pure Neumann not yet implemented') else: solver = pg.core.LinSolver(False) solver.setMatrix(A, 0) if singleForce: if isComplex is True: # clean this up rhs = pg.core.toComplex(rhs.real, rhs.imag) u = solver.solve(rhs).array() else: u = solver.solve(rhs) else: for i, r in enumerate(rhs): u[i] = solver.solve(r) solverTime = swatch.duration(True) if verbose: if stats: stats.solverTime = solverTime print("Solving time: ", solverTime) if len(kwargs.keys()) > 0: pg.warn("Unused arguments", *kwargs) return u else: # times given pg.solver.checkCFL(times, mesh, max(np.array(a).flatten())) if debug: print("start TL", swatch.duration()) if c != 1.0: c = cellValues(mesh, c, userData=userData) M = createMassMatrix(mesh, c) F = createLoadVector(mesh, f) u0 = np.zeros(dof) if 'u0' in kwargs: u0 = parseArgToArray(kwargs['u0'], dof, mesh, userData) progress = None if 'progress' in kwargs: from pygimli.utils import ProgressBar progress = ProgressBar(its=len(times), width=40, sign='+') theta = kwargs.pop('theta', 1.0) dynamic = kwargs.pop('dynamic', False) if not dynamic: S = createStiffnessMatrix(mesh, a) assembleBC(bc, mesh, S, F, time=0.0, userData=userData) return crankNicolson(times, S, M, f=F, u0=u0, theta=theta, progress=progress) rhs = np.zeros((len(times), dof)) # no time dependency for rhs so far ... TODO rhs[:] = F # this is slow: optimize if debug: print("rhs", swatch.duration()) U = np.zeros((len(times), dof)) U[0, :] = u0 # init state u = pg.Vector(dof, 0.0) if debug: print("u0", swatch.duration()) measure = 0. for n in range(1, len(times)): swatch.reset() dt = times[n] - times[n-1] # previous timestep # print "i: ", i, dt, U[i - 1] swatch.reset() # (A + a*B)u is fastest, # followed by A*u + (B*u)*a and finally A*u + a*B*u and br = (M + S*(dt * (theta - 1.))) * U[n - 1] + \ dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n]) # print ('a',swatch.duration(True)) # br = M * U[n - 1] - (A * U[n - 1]) * (dt*(1.0 - theta)) + \ # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n]) # print ('br',swatch.duration(True)) # br = M * U[n - 1] - (dt*(1.0 - theta)) * A * U[n - 1] + \ # dt * ((1.0 - theta) * rhs[n - 1] + theta * rhs[n]) # print ('c',swatch.duration(True)) measure += swatch.duration() A = M + S * dt * theta assembleBC(bc, mesh, A, br, time=times[n], userData=userData) if 'assembleOnly' in kwargs: return A, br # u = S/b t_prep = swatch.duration(True) solver = pg.core.LinSolver(A, verbose) solver.solve(br, u) if 'plotTimeStep' in kwargs: kwargs['plotTimeStep'](u, times[n]) U[n, :] = np.asarray(u) if progress: progress.update(n, 't_prep: {0}ms t_step {1}s'.format( pg.pf(t_prep*1000), pg.pf(swatch.duration()))) if debug: print("Measure(" + str(len(times)) + "): ", measure, measure / len(times)) return U
[docs] def checkCFL(times, mesh, vMax, verbose=False): """Check Courant-Friedrichs-Lewy condition. For advection and flow problems. CFL Number should be lower then 1 to ensure stability. Parameters ---------- """ if pg.isScalar(times): dt = times else: dt = times[1] - times[0] dx = min(mesh.h()) # min(entity.shape().h() # if mesh.dimension() == 1: # dx = min(mesh.cellSizes()) # else: # dx = min(mesh.boundarySizes()) c = vMax * dt / dx if c > 1: pg.warn("Courant-Friedrichs-Lewy Number:", c, "but should be lower 1 to ensure movement inside a cell " "per timestep. (" "vMax =", vMax, "dt =", dt, "dx =", dx, "dt <", dx/vMax, " | N > ", int(dt/(dx/vMax))+1, ")") if verbose: pg.info("Courant-Friedrichs-Lewy Number:", c) return c
[docs] def crankNicolson(times, S, I, f=None, u0=None, theta=1.0, dirichlet=None, solver=None, progress=None): """Generic Crank Nicolson solver for time dependend problems. Limitations so far: S = Needs to be constant over time (i.e. no change in coefficients) f = constant over time (would need assembling in every step) Args ---- times: iterable(float) Timeteps to solve for. Give at least 2. S: Matrix Systemmatrix holds your discrete equations and boundary conditions I: Matrix Identity matrix (FD, FV) or Masselementmatrix (FE) to handle solution vector u0: iterable [None] Starting condition. zero if not given f: iterable (float) [None] External forces. Note f might also contain compensation values due to algebraic Dirichlet correction of S theta: float [1.0] * 0: Backward difference scheme (implicit) * 1: Forward difference scheme (explicit) strong time steps dependency .. will be unstable for to small values * 0.5: probably best tradeoff but can also be unstable dirichlet: dirichlet generator Genertor object to applay dirichlet boundary conditions solver: LinSolver [None] Provide a pre configured solver if you want some special. progress: Progress [None] Provide progress object if you want to see some. Returns ------- np.ndarray: Solution for each time steps """ if len(times) < 2: raise BaseException("We need at least 2 times for " "Crank-Nicolsen time discretization." + str(len(times))) # sw = pg.core.Stopwatch(True) timeAssemble = [] timeSolve = [] timeMeasure = False if progress: timeMeasure = True dof = S.rows() rhs = np.zeros((len(times), dof)) if f is not None: rhs[:] = f u = np.zeros((len(times), dof)) if u0 is not None: u[0, :] = u0 if theta == 0: A = I.copy() if solver is None: solver = pg.solver.LinSolver(solver='scipy') dt = 0.0 for n in range(1, len(times)): newDt = times[n] - times[n-1] if abs(newDt - dt) > 1e-8: # new dt, so we need to factorize the matrix again dt = newDt # pg.info('dt', dt) A = I + S * (dt * theta) if dirichlet is not None: dirichlet.apply(A) solver.factorize(A) St = None if timeMeasure: pg.tic(key='CrankNicolsonLoop') if theta == 0: if St is None: St = I - S * dt # cache what's possible b = St * u[n-1] + dt * rhs[n-1] elif theta == 1: b = I * u[n-1] + dt * rhs[n] else: if St is None: St = I - S * (dt*(1.-theta)) # cache what's possible b = St * u[n-1] + dt * ((1.0 - theta) * rhs[n-1] + theta * rhs[n]) if dirichlet is not None: dirichlet.apply(b) if timeMeasure: timeAssemble.append(pg.dur(key='CrankNicolsonLoop', reset=True)) u[n, :] = solver(b) if timeMeasure: timeSolve.append(pg.dur(key='CrankNicolsonLoop')) if progress: progress.update( n, 't_prep: ' + pg.pf(timeAssemble[-1]*1000) + 'ms ' + 't_step: ' + pg.pf(timeSolve[-1]*1000) + 'ms') # if verbose and (n % verbose == 0): # # print(min(u[n]), max(u[n])) # print("timesteps:", n, "/", len(times), # 'runtime:', sw.duration(), "s", # 'assemble:', np.mean(timeAssemble), # 'solve:', np.mean(timeSolve)) return u
[docs] class RungeKutta(object): """TODO DOCUMENT ME""" rk4a = [0.0, -567301805773.0/1357537059087.0, -2404267990393.0/2016746695238.0, -3550918686646.0/2091501179385.0, -1275806237668.0/842570457699.0] rk4b = [1432997174477.0/9575080441755.0, 5161836677717.0/13612068292357.0, 1720146321549.0/2090206949498.0, 3134564353537.0/4481467310338.0, 2277821191437.0/14882151754819.0] rk4c = [0.0, 1432997174477.0/9575080441755.0, 2526269341429.0/6820363962896.0, 2006345519317.0/3224310063776.0, 2802321613138.0/2924317926251.0]
[docs] def __init__(self, solver, verbose=False): """TODO DOCUMENT_ME""" self.solver = solver self.verbose = verbose self.order = 5 self.dt = None self.time = 0 self.tMax = None self.u = None self.resU = None self.nSteps = 0
[docs] def run(self, u0, dt, tMax=1): """TODO DOCUMENT_ME""" self.start(u0, dt, tMax) for _ in range(self.nSteps): self.step() return self.u
[docs] def start(self, u0, dt, tMax=1): """TODO DOCUMENT_ME""" self.nSteps = int(np.ceil(tMax/dt)) self.dt = dt self.time = 0 self.tMax = tMax self.u = deepcopy(u0) self.resU = deepcopy(u0) if isinstance(self.resU, list): for r in self.resU: r *= 0.0 else: self.resU *= 0.0
[docs] def step(self): """TODO DOCUMENT ME""" if self.time + self.dt > self.tMax: self.dt = self.tMax - self.time if self.order == 1: # explicit Euler k1 = self.solver.explicitRHS(self.u, self.time) self.u += self.dt * k1 elif self.order == 3: k1 = self.solver.explicitRHS(self.u, self.time) k1 = self.u + self.dt * k1 k2 = self.solver.explicitRHS(k1, self.time) k2 = (3*self.u + k1 + self.dt*k2)/4 k3 = self.solver.explicitRHS(k2, self.time) self.u = (self.u + 2*k2 + 2*self.dt*k3)/3 elif self.order == 4: # classical 4 step Runge-Kutta rk4 k1 = self.solver.explicitRHS(self.u, self.time) k2 = self.solver.explicitRHS(self.u + self.dt/2 * k1, self.time + self.dt/2) k3 = self.solver.explicitRHS(self.u + self.dt/2 * k2, self.time + self.dt/2) k4 = self.solver.explicitRHS(self.u + self.dt * k3, self.time + self.dt) self.u += 1./6. * self.dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) elif self.order == 5: # low storage Version of rk4 for jRK in range(5): tLocal = self.time + self.rk4c[jRK] * self.dt rhs = self.solver.explicitRHS(self.u, tLocal) if isinstance(self.resU, list): for i in range(len(self.resU)): self.resU[i] = self.rk4a[jRK] * self.resU[i] + \ self.dt * rhs[i] self.u[i] += self.rk4b[jRK] * self.resU[i] else: self.resU = self.rk4a[jRK] * self.resU + self.dt * rhs self.u += self.rk4b[jRK] * self.resU self.time += self.dt return self.u
if __name__ == "__main__": pass