pygimli.solver

General physics independent solver interface.

Overview

Functions

anisotropyMatrix(lon, trans, theta)

Create anisotropy matrix with desired properties.

applyBoundaryValues(uB, mesh, uBBC)

TODO Documentme.

assembleBC(bc, mesh, mat, rhs[, time, …])

Shortcut to apply all boundary conditions.

assembleDirichletBC(mat, boundaryPairs[, …])

Apply Dirichlet boundary condition.

assembleLoadVector(mesh, f[, userData])

Assemble the load vector.

assembleNeumannBC(rhs, boundaryPairs[, …])

Apply Neumann condition to the system matrix S.

assembleRobinBC(mat, boundaryPairs[, rhs, …])

Apply Robin boundary condition.

boundaryNormals(mesh)

Collect all boundary outer normal vectors.

cellDataToBoundaryData(mesh, v)

TODO Documentme.

cellDataToBoundaryDataMatrix(mesh)

TODO Documentme.

cellDataToBoundaryGrad(mesh, v, vGrad)

TODO Documentme.

cellDataToCellGrad(mesh, v, CtB)

TODO Documentme.

cellDataToCellGrad2(mesh, v)

TODO Documentme.

cellToFace(boundary, vec[, harmonic])

DEPRECATED.

cellToFaceArithmetic(boundary, AMM)

TODO Documentme.

cellValues(mesh, arg, **kwargs)

Get a value for each cell.

checkCFL(times, mesh, vMax)

Check Courant-Friedrichs-Lewy condition.

constitutiveMatrix([lam, mu, E, nu, dim, …])

Create constitutive matrix for 2 or 3D isotropic media.

crankNicolson(times, theta, matS, matI, f[, …])

S = constant over time f = constant over time

createFVPostProzessMesh(mesh, u, uDirichlet)

Create node based post processing of cell centered FV solutions.

createForceVector(mesh, f[, userData])

Create a right hand side vector for vector valued solutions.

createLoadVector(mesh[, f, userData])

Create right hand side vector based on the given mesh and load values (scalar solution) or force vectors (vector value solution).

createMassMatrix(mesh[, b])

Create the mass matrix.

createMesh(poly[, quality, area, smooth, …])

Create a mesh for a given PLC or point list.

createStiffnessMatrix(mesh[, a, isVector])

Create the Stiffness matrix.

deepcopy(x[, memo, _nil])

Deep copy operation on arbitrary Python objects.

diffusionConvectionKernel(mesh[, a, b, uB, …])

Generate system matrix for diffusion and convection in a velocity field.

div(mesh, v)

Return the discrete interpolated divergence field.

divergence(mesh[, func, normMap, order])

Divergence for callable function func((x,y,z)).

findDiffusion(mesh, a, b, c[, nc])

TODO Documentme.

findVelocity(mesh, v, b, c[, nc])

Find velocity for boundary b based on vector field v.

generateBoundaryValue(boundary, arg[, time, …])

Generate a value for the given Boundary.

grad(mesh, u[, r])

Return the discrete interpolated gradient \(\mathbf{v}\) for a given scalar field \(\mathbf{u}\).

greenDiffusion1D(x[, t, a, dim])

Greens function for diffusion operator.

identity(dom[, start, end, scale])

Create identity matrix.

intDomain(u[, mesh])

Return integral over nodal solution \(u\).

linSolve(mat, b[, solver, verbose])

Direct solution after \(\textbf{x}\) using core LinSolver.

normH1(u[, mat, mesh])

Create (H1) norm for the finite element space.

normL2(u[, mat, mesh])

Create Lebesgue (L2) norm for finite element space.

parseArgPairToBoundaryArray(pair, mesh)

Parse boundary related pair argument to create a list of [ GIMLI::Boundary, value|callable ].

parseArgToArray(arg, nDof[, mesh, userData])

Parse array related arguments to create a valid value array.

parseArgToBoundaries(args, mesh)

Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ]

parseDictKey_(key, markers)

Parse dictionary key of type str to marker list.

parseMapToCellArray(attributeMap, mesh[, …])

Parse a value map to cell attributes.

show([obj, data])

Mesh and model visualization.

showSparseMatrix(mat[, full])

Show the content of a sparse matrix.

solve(mesh, **kwargs)

Solve partial differential equation.

solveFiniteElements(mesh[, a, b, f, bc, …])

Solve partial differential equation with Finite Elements.

solveFiniteVolume(mesh[, a, b, f, fn, vel, …])

Solve partial differential equation with Finite Volumes.

triDiagToeplitz(dom, a, l, r[, start, end])

Create tri-diagonal Toeplitz matrix.

Classes

RungeKutta(solver[, verbose])

TODO DOCUMENT ME

WorkSpace

Functions

anisotropyMatrix

pygimli.solver.anisotropyMatrix(lon, trans, theta)[source]

Create anisotropy matrix with desired properties.

Anistropy tensor from longitudinal value lon, transverse value trans and the angle theta of the symmetry axis relative to the vertical after cite:WieseGreZho2015 https://www.researchgate.net/publication/249866312_Explicit_expressions_for_the_Frechet_derivatives_in_3D_anisotropic_resistivity_inversion

Examples using pygimli.solver.anisotropyMatrix

applyBoundaryValues

pygimli.solver.applyBoundaryValues(uB, mesh, uBBC)[source]

TODO Documentme.

assembleBC

pygimli.solver.assembleBC(bc, mesh, mat, rhs, time=None, userData={}, dofOffset=0)[source]

Shortcut to apply all boundary conditions.

Shortcut to apply all boundary conditions will only forward to appropriate assemble functions.

Returns

Return type

None

assembleDirichletBC

pygimli.solver.assembleDirichletBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, nodePairs=None, dofOffset=0)[source]

Apply Dirichlet boundary condition.

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.

\[u(\textbf{r}, t) = h \quad\text{for}\quad\textbf{r}\quad\text{on}\quad\delta\Omega=\Gamma_{\text{Dirichlet}}\]
Parameters
  • mat (GIMLI::RSparseMatrix) – System matrix of the system equation.

  • boundaryPair (list()) –

    List of pairs [GIMLI::Boundary, h]. The value \(h\) will assigned to the nodes of the boundaries. Later assignment overwrites prior.

    \(h\) need to be a scalar value (float or int) or a value generator function that will be executed at runtime. See pygimli.solver.solver.parseArgToBoundaries and Modelling with Boundary Conditions 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.

  • rhs (GIMLI::RVector) – Right hand side vector of the system equation will bet set to \(u_{\text{D}}\)

  • time (float) – Will be forwarded to value generator.

  • userData (class) – Will be forwarded to value generator.

  • dofOffset (int[0]) – Offset for matrix index.

Examples using pygimli.solver.assembleDirichletBC

assembleLoadVector

pygimli.solver.assembleLoadVector(mesh, f, userData={})[source]

Assemble the load vector. See createLoadVector. Maybe we will remove this

assembleNeumannBC

pygimli.solver.assembleNeumannBC(rhs, boundaryPairs, nDim=1, time=0.0, userData={}, dofOffset=0)[source]

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.

\[\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 (Vector) – Right hand side vector of length node count.

  • boundaryPair (list()) –

    List of pairs [ GIMLI::Boundary, g ]. The value \(g\) will assigned to the nodes of the boundaries. Later assignment overwrites prior.

    \(g\) need to be a scalar value (float or int) or a value generator function that will be executed at run time.

    See pygimli.solver.solver.parseArgToBoundaries and Modelling with Boundary Conditions for example syntax,

  • nDim (int [1]) – Number of dimensions for vector valued problems. The rhs array need top 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.

assembleRobinBC

pygimli.solver.assembleRobinBC(mat, boundaryPairs, rhs=None, time=0.0, userData={}, dofOffset=0)[source]

Apply Robin boundary condition.

Apply Robin boundary condition to the system matrix and the rhs vector

\[\begin{split}\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}}\\\end{split}\]
Parameters
  • mat (GIMLI::SparseMatrix) – System matrix of the system equation.

  • boundaryPair (list()) –

    List of pairs [GIMLI::Boundary, \(a, u_0\) | \(\alpha, \beta, \gamma\)]. The values will assigned to the nodes of the boundaries. Later assignment overwrites prior.

    Values can be a single value for \(\alpha\) or \(a\), two values will be interpreted as \(a, u_0\), and three values will be \(\alpha, \beta, \gamma\). Also generator (callable) is possible which will be executed at run time. See pygimli.solver.solver.parseArgToBoundaries Modelling with Boundary Conditions 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.

boundaryNormals

pygimli.solver.boundaryNormals(mesh)[source]

Collect all boundary outer normal vectors.

cellDataToBoundaryData

pygimli.solver.cellDataToBoundaryData(mesh, v)[source]

TODO Documentme.

cellDataToBoundaryDataMatrix

pygimli.solver.cellDataToBoundaryDataMatrix(mesh)[source]

TODO Documentme.

cellDataToBoundaryGrad

pygimli.solver.cellDataToBoundaryGrad(mesh, v, vGrad)[source]

TODO Documentme.

cellDataToCellGrad

pygimli.solver.cellDataToCellGrad(mesh, v, CtB)[source]

TODO Documentme.

cellDataToCellGrad2

pygimli.solver.cellDataToCellGrad2(mesh, v)[source]

TODO Documentme.

cellToFace

pygimli.solver.cellToFace(boundary, vec, harmonic=False)[source]

DEPRECATED.

Interpolate cell to face values by weighted arithmetic/harmonic mean.

cellToFaceArithmetic

pygimli.solver.cellToFaceArithmetic(boundary, AMM)[source]

TODO Documentme.

cellValues

pygimli.solver.cellValues(mesh, arg, **kwargs)[source]

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

Parameters
  • mesh (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, the dict 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 parseDictKey.

    Iterable of length mesh.nodeCount() will be interpolated to cellCenters.

  • userData (class) – Used if arg is callable

Returns

ret – Array of desired length filled with the appropriate values.

Return type

GIMLI::RVector | ndarray(mesh.cellCount(), xx )

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.]]

 [[1. 0.]
  [0. 2.]]

 [[1. 0.]
  [0. 2.]]

 [[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]]

 [[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]

 [[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]

 [[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]

checkCFL

pygimli.solver.checkCFL(times, mesh, vMax)[source]

Check Courant-Friedrichs-Lewy condition.

For advection and flow problems. CFL Number should be lower then 1 to ensure stability.

constitutiveMatrix

pygimli.solver.constitutiveMatrix(lam=None, mu=None, E=None, nu=None, dim=2, voigtNotation=False)[source]

Create constitutive matrix for 2 or 3D isotropic media.

Either give lam and mu or E and nu.

Parameters
  • lam (float [None]) –

    1. Lame’ constant

  • mu (float [None]) –

    1. Lame’ constant = G = Schubmodul

  • E (float [None]) – Young’s Modulus

  • nu (float [None]) – Poisson’s ratio

  • voigtNotation (bool [False]) – Return in Voigt’s notation instead of Kelvin’s notation [default].

Returns

C – Either 3x3 or 6x6 matrix depending on the dimension

Return type

mat

crankNicolson

pygimli.solver.crankNicolson(times, theta, matS, matI, f, u0=None, progress=None)[source]

S = constant over time f = constant over time

createFVPostProzessMesh

pygimli.solver.createFVPostProzessMesh(mesh, u, uDirichlet)[source]

Create node based post processing of cell centered FV solutions.

Create a mesh suitable for node based post processing of cell centered Finite Volume solutions. This is something like cellDataToNodeData with extra Dirichlet points but without smoothing due to interpolation.

IMPROVE DOC!!

createForceVector

pygimli.solver.createForceVector(mesh, f, userData={})[source]

Create a right hand side vector for vector valued solutions.

Parameters
  • f ([ convertable ]) – List of rhs side options. Must be convertable to createLoadVector. See createLoadVector

  • rhs (np.array()) – Squeezed vector of length mesh.nodeCount() * mesh.dimensions()

createLoadVector

pygimli.solver.createLoadVector(mesh, f=1.0, userData={})[source]

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.

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() will be assumed to be allready 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 – Right hand side load vector for scalar values or squeezed vector values.

Return type

pg.Vector(mesh.nodeCount())

Examples using pygimli.solver.createLoadVector

createMassMatrix

pygimli.solver.createMassMatrix(mesh, b=None)[source]

Create the mass matrix.

Calculates the Mass matrix (Finite element identity matrix) the given mesh.

..math::

Parameters
  • mesh (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 – Mass element matrix

Return type

GIMLI::RSparseMatrix

createMesh

pygimli.solver.createMesh(poly, quality=32, area=0.0, smooth=None, switches=None, verbose=False, **kwargs)[source]

Create a mesh for a given PLC or point list.

The mesh is created by triangle or tetgen if the pgGIMLi support for these mesh generators is installed. A PLC needs to contain nodes and boundaries and should be valid in the sense that the boundaries are non-intersecting.

If poly is a list of coordinates, a simple Delaunay mesh with a convex hull will be created. Quality and area arguments are ignored for this case to create a mesh with one node for each coordinate position.

Parameters
  • poly (GIMLI::Mesh or list or ndarray) –

    • 2D or 3D gimli mesh that contains the PLC.

    • 2D mesh needs edges

    • 3D mesh needs a plc and tetgen as system component

    • List of x y pairs [[x0, y0], … ,[xN, yN]]

    • ndarray [x_i, y_i]

    • PLC or list of PLCs

  • quality (float) – 2D triangle quality sets a minimum angle constraint. Be careful with values above 34 degrees. 3D tetgen quality. Be careful with values below 1.12.

  • area (float) – Maximum element size (global). 2D maximum triangle size in m*², 3D maximum tetrahedral size in m³.

  • smooth (tuple) – [smoothing algorithm, number of iterations] 0: no smoothing 1: node center 2: weighted node center

  • switches (str) – Set additional triangle command switches. https://www.cs.cmu.edu/~quake/triangle.switch.html

Returns

mesh

Return type

GIMLI::Mesh

>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>> rect = mt.createRectangle(start=[0, 0], end=[4, 1])
>>> ax, _ = pg.show(mt.createMesh(rect, quality=10))
>>> ax, _ = pg.show(mt.createMesh(rect, quality=33))
>>> ax, _ = pg.show(mt.createMesh(rect, quality=33, area=0.01))
>>> pg.wait()
../../_images/pygimli-solver-1_00.png

(png, pdf)

../../_images/pygimli-solver-1_01.png

(png, pdf)

../../_images/pygimli-solver-1_02.png

(png, pdf)

createStiffnessMatrix

pygimli.solver.createStiffnessMatrix(mesh, a=None, isVector=False)[source]

Create the Stiffness matrix.

Calculates the Stiffness matrix \({\bf S}\) for the given mesh scaled with the per cell values a.

..math::

Parameters
  • mesh (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 – Stiffness matrix, with real or complex values.

Return type

GIMLI::[C]SparseMatrix | [C]SparseMapMatrix

Examples using pygimli.solver.createStiffnessMatrix

deepcopy

pygimli.solver.deepcopy(x, memo=None, _nil=[])[source]

Deep copy operation on arbitrary Python objects.

See the module’s __doc__ string for more info.

diffusionConvectionKernel

pygimli.solver.diffusionConvectionKernel(mesh, a=None, b=0.0, uB=None, duB=None, vel=0, fn=None, scheme='CDS', sparse=False, time=0.0, userData=None)[source]

Generate system matrix for diffusion and convection in a velocity field.

Particle concentration u inside a velocity field.

Peclet Number - ratio between convection/diffusion = F/D

F = velocity flow trough volume boundary, D = diffusion coefficient

Parameters
  • mesh (GIMLI::Mesh) – Mesh represents spatial discretization of the calculation domain

  • a (value | array | callable(cell, userData)) – Diffusion coefficient per cell

  • b (value | array | callable(cell, userData)) – TODO What is b

  • fn (iterable(cell)) – TODO What is fn

  • vel (ndarray (N,dim) | RMatrix(N,dim)) – velocity field [[v_i,]_j,] with i=[1..3] for the mesh dimension and j = [0 .. N-1] per Cell or per Node so N is either mesh.cellCount() or mesh.nodeCount()

  • scheme (str [CDS]) –

    Finite volume scheme * CDS – Central Difference Scheme.

    maybe irregular for Peclet no. |F/D| > 2 Diffusion dominant. Error of order 2

    • UDS – Upwind Scheme.

      Convection dominant. Error of order 1

    • HS – Hybrid Scheme.

      Diffusion dominant for Peclet-number |(F/D)| < 2 Convection dominant else.

    • PS – Power Law Scheme.

      Identical to HS for Peclet-number |(F/D)| > 10 and near to ES else Convection dominant.

    • ES – Exponential scheme

      Only stationary one-dimensional but exact solution

Returns

  • S (GIMLI::SparseMatrix | numpy.ndarray(nCells, nCells)) – Kernel matrix, depends on vel, a, b, scheme, uB, duB .. if some of this has been changed you cannot cache these matrix

  • rhsBoundaryScales (ndarray(nCells)) – RHS offset vector

div

pygimli.solver.div(mesh, v)[source]

Return the discrete interpolated divergence field.

Return the discrete interpolated divergence field. \(\mathbf{u}\) for each cell for a given vector field \(\mathbf{v}\). First order integration via boundary center.

\[\begin{split}d(cells) & = \nabla\cdot\vec{v} \\ d(c_i) & = \sum_{j=0}^{N_B}\vec{v}_{B_j} \cdot \vec{n}_{B_j}\end{split}\]
Parameters
  • mesh (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 of divergence values for each cell in the given mesh.

Return type

array(M)

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 are wrong due to interpolation to the 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

divergence

pygimli.solver.divergence(mesh, func=None, normMap=None, order=1)[source]

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.

findDiffusion

pygimli.solver.findDiffusion(mesh, a, b, c, nc=None)[source]

TODO Documentme.

Parameters
  • a – Attribute for diffusion coefficient. Cell based or Boundary based

  • b (GIMLI::Boundary) – Boundary

  • c (GIMLI::Cell) – Associated cell in flow direction

  • nc (GIMLI::Cell) – associated neighbor cell .. if one given from flow direction

findVelocity

pygimli.solver.findVelocity(mesh, v, b, c, nc=None)[source]

Find velocity for boundary b based on vector field v.

Parameters
  • mesh (GIMLI::Mesh) –

  • v (array [(N,dim)]) – velocity field [[v_i,]_j,] with i=[1..3] for the mesh dimension and j = [0 .. N-1] per Cell or per Node so N is either mesh.cellCount() or mesh.nodeCount()

  • b (GIMLI::Boundary) – Boundary

  • c (GIMLI::Cell) – Associated Cell in flow direction

  • nc (GIMLI::Cell) – Associated neighbor cell .. if one given from flow direction

generateBoundaryValue

pygimli.solver.generateBoundaryValue(boundary, arg, time=0.0, userData={}, expectList=False)[source]

Generate a value for the given Boundary.

Parameters
  • boundary (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 need to return a ndarray((nodes, dims)) like ‘lambda _b: np.array([[force_x, force_y, force_z] for n in _b.nodes()]).T’

Returns

val – Value for all nodes of the boundary.

Return type

[float]

grad

pygimli.solver.grad(mesh, u, r=None)[source]

Return the discrete interpolated gradient \(\mathbf{v}\) for a given scalar field \(\mathbf{u}\).

\[\begin{split}\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}}\end{split}\]

With \(\mathcal{N}=\cup_{i=0}^{N} \text{Node}_i\), \(\mathcal{C}=\cup_{j=0}^{M} \text{Cell}_j\), \(\mathbf{u}=\{u(\mathbf{r}_{i})\}\in I\!R\) and \(\mathbf{r}_{i} = (x_i, y_i, z_i)^{\text{T}}\)

The discrete scalar field \(\mathbf{u}(\mathbf{r}_{\mathcal{N}})\) need to be defined for each node position \(\mathbf{r}_{\mathcal{N}}\). The resulting vector field \(\mathbf{v}(\mathbf{r}_{\mathcal{C}})\) is defined for each cell center position \(\mathbf{r}_{\mathcal{C}}\). If you need other positions than the cell center, provide an appropriate array of coordinates \(\mathbf{r}\).

Parameters
  • mesh (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 – Resulting vector field defined on \(\mathbf{v}(\mathbf{r}_{\mathcal{C}})\). M is number of cells or length of given alternative coordinates r.

Return type

ndarray((M, 3))

>>> 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()

(png, pdf)

../../_images/pygimli-solver-2.png

Examples using pygimli.solver.grad

greenDiffusion1D

pygimli.solver.greenDiffusion1D(x, t=0, a=1, dim=1)[source]

Greens function for diffusion operator.

Provides the elementary solution for:

\[g(x,t) = \partial t + a \Delta\]

To find a solution for:

\[u(x,t) \quad\text{for}\quad\frac{\partial u(x,t)}{\partial t} + a \Delta u(x,t) = f(x)\]
\[x = [-x, 0, x]\]
\[u(x,t) = g(x,t) * f(x)\]
Parameters
  • x (array_like) – Discrete spatial coordinates. Should better be symmetric [-x, 0, x].

  • t (float) – Time

  • a (float, optional) – Scale for Laplacian operator

  • dim (int, optional) – Spatial dimension [1]

Returns

g – Discrete Greens’function

Return type

array_like

>>> import numpy as np
>>> from pygimli.solver import greenDiffusion1D
>>> dx = 0.001
>>> x = np.arange(0, 0.1+dx, dx)
>>> g = greenDiffusion1D(np.hstack((-x[:0:-1], x)), t=1.0, a=1e-4)
>>> g *= dx
>>> f = np.zeros(len(x))
>>> f[int(len(x)/2)] = 1.
>>> u = np.convolve(g, f)[len(x)-1:2*len(x)-1]
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> _ = ax.plot(x, u, label="u(x)")
>>> _ = ax.plot(x, g[::2], label="g(x)")
>>> _ = ax.legend()
>>> fig.show()

(png, pdf)

../../_images/pygimli-solver-3.png

identity

pygimli.solver.identity(dom, start=0, end=- 1, scale=1)[source]

Create identity matrix.

intDomain

pygimli.solver.intDomain(u, mesh=None)[source]

Return integral over nodal solution \(u\).

\[\int_{\Omega} u\]

linSolve

pygimli.solver.linSolve(mat, b, solver=None, verbose=False)[source]

Direct solution after \(\textbf{x}\) using core LinSolver.

\[\textbf{A}\textbf{x} = \textbf{b}\]

If \(\textbf{A}\) is symmetric, sparse and positive definite.

Parameters
  • mat (GIMLI::RSparseMatrix | 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 – Solution vector

Return type

GIMLI::RVector

Examples using pygimli.solver.linSolve

normH1

pygimli.solver.normH1(u, mat=None, mesh=None)[source]

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 (GIMLI::Mesh) – Mesh with the FE space to generate S if necessary.

Returns

ret\(H1(u)\) norm.

Return type

float

normL2

pygimli.solver.normL2(u, mat=None, mesh=None)[source]

Create Lebesgue (L2) norm for finite element space.

Find the L2 Norm for a solution for the finite element space. \(u\) exact solution \({\bf M}\) Mass matrix, i.e., Finite element identity matrix.

\[\begin{split}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}}\end{split}\]

The error for any approximated solution \(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 (GIMLI::Mesh) – Mesh with the FE space to generate M if necessary.

Returns

ret\(L2(u)\) norm.

Return type

float

Examples using pygimli.solver.normL2

parseArgPairToBoundaryArray

pygimli.solver.parseArgPairToBoundaryArray(pair, mesh)[source]

Parse boundary related pair argument to create a list of [ 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 pygimli.solver.solver.generateBoundaryValue and distributed to each boundary. Callable functions will be executed at run time. ‘*’ will be interpreted as all boundary elements with one neighboring cell

  • mesh (GIMLI::Mesh) – Used to find boundaries by marker.

Returns

bc – [GIMLI::Boundary, value|callable]

Return type

list()

parseArgToArray

pygimli.solver.parseArgToArray(arg, nDof, mesh=None, userData={})[source]

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 GIMLI::Cell , GIMLI::Node or GIMLI::Boundary depending on nDof, where nDof is mesh.cellCount(), mesh.nodeCount() or mesh.boundaryCount(), respectively.

  • nDof (int | [int]) – Desired array size.

  • mesh (GIMLI::Mesh) – Used if arg is callable

  • userData (class) – Used if arg is callable

Returns

ret – Array of desired length filled with the appropriate values.

Return type

GIMLI::RVector

Examples using pygimli.solver.parseArgToArray

parseArgToBoundaries

pygimli.solver.parseArgToBoundaries(args, mesh)[source]

Parse boundary related arguments to create a valid boundary value list: [ GIMLI::Boundary, value|callable ]

Parameters
  • args (dict, float, callable) –

    Dictionary is preferred (key=value|callable). If args is just a callable or float every outer boundary will be 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 pygimli.solver.solver.parseArgPairToBoundaryArray

  • mesh (GIMLI::Mesh) – Used to find boundaries by marker

Returns

boundaries – [ GIMLI::Boundary, value|callable ]

Return type

list()

>>> # 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()

(png, pdf)

../../_images/pygimli-solver-4.png

Examples using pygimli.solver.parseArgToBoundaries

parseDictKey

pygimli.solver.parseDictKey_(key, markers)[source]

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) – Supported are - ‘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 – List of integers described by key

Return type

[int]

parseMapToCellArray

pygimli.solver.parseMapToCellArray(attributeMap, mesh, default=0.0)[source]

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 (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 of length mesh.cellCount()

Return type

array

Examples using pygimli.solver.parseMapToCellArray

show

pygimli.solver.show(obj=None, data=None, **kwargs)[source]

Mesh and model visualization.

Syntactic sugar to show a obj with data. Forwards to a known visualization for obj. Typical is pygimli.viewer.showMesh or pygimli.viewer.mayaview.showMesh3D to show most of the typical 2D and 3D content. See tutorials and examples for usage hints. An empty show call creates an empty ax window.

Parameters
  • obj (obj) – obj can be so far. * GIMLI::Mesh or list of meshes * DataContainer * pg.core.Sparse[Map]Matrix

  • data (iterable) – Optionally data to visualize. See appropriate show function.

Keyword Arguments

**kwargs

Additional kwargs forward to appropriate show functions.

  • axaxe [None]

    Matplotlib axes object. Create a new if necessary.

  • fitViewbool [True]

    Scale x and y limits to match the view.

Returns

  • Return the results from the showMesh functions. Usually the axe object*

  • and a colorbar.

See also

showMesh()

showSparseMatrix

pygimli.solver.showSparseMatrix(mat, full=False)[source]

Show the content of a sparse matrix.

Parameters

solve

pygimli.solver.solve(mesh, **kwargs)[source]

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 pygimli.solver.solveFiniteElements

TODO pygimli.solver.solveFiniteVolume

solveFiniteElements

pygimli.solver.solveFiniteElements(mesh, a=1.0, b=None, f=0.0, bc=None, times=None, userData={}, verbose=False, **kwargs)[source]

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:

\[\begin{split}\frac{\partial u}{\partial t} & = \nabla\cdot(a \nabla u) + b u + f(\mathbf{r},t)~~|~~\Omega_{\text{Mesh}}\\ u & = h~~|~~\Gamma_{\text{Dirichlet}}\\ \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}}\end{split}\]

for the scalar \(u(\mathbf{r}, t)\) or vector \(\mathbf(u)(\mathbf{r}, t)\) solution at each node of a given mesh. The Domain \(\Omega\) and the Boundary \(\Gamma\) are defined through the mesh with appropriate boundary marker.

Note, to ensure vector solution either set vector forces or at least on vector component boundary condition.

Parameters
  • mesh (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, anisotropic matrix or elastic tensor.

  • b (value | array | callable(cell, userData) [None]) – Cell values. None means the term is unused.

  • 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 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]}.

  • 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]
  • \(theta = 0\) means explicit Euler, maybe stable for

\(\Delta t \quad\text{near}\quad h\) - \(theta = 0.5\), Crank-Nicolson, maybe instable - \(theta = 1\), implicit Euler

If unsure choose \(\theta = 0.5 + \epsilon\), which is 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.

pureNeumann: bool [auto]

If set or detected automatic, we add the additional condition: \(\int_domain u dv = 0\) which makes elliptic problems well posed again.

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 will be used for the resulting array.

Returns

u – Returns the solution u either 1,n array for stationary problems or for m,n array for m time steps

Return type

array

>>> import pygimli as pg
>>> from pygimli.meshtools import polytools as plc
>>> from pygimli.viewer.mpl import drawField, drawMesh
>>> import matplotlib.pyplot as plt
>>> world = plc.createWorld(start=[-10, 0], end=[10, -10],
...                         marker=1, worldMarker=False)
>>> c1 = plc.createCircle(pos=[0.0, -5.0], radius=3.0, area=.1, marker=2)
>>> mesh = pg.meshtools.createMesh([world, c1], quality=34.3)
>>> u = pg.solver.solveFiniteElements(mesh, a={1: 100.0, 2: 1.0},
...                                   bc={'Dirichlet':{4: 1.0, 2: 0.0}})
>>> fig, ax = plt.subplots()
>>> pc = drawField(ax, mesh, u)
>>> drawMesh(ax, mesh)
>>> plt.show()

(png, pdf)

../../_images/pygimli-solver-5.png

Examples using pygimli.solver.solveFiniteElements

solveFiniteVolume

pygimli.solver.solveFiniteVolume(mesh, a=1.0, b=0.0, f=0.0, fn=0.0, vel=None, u0=0.0, times=None, uL=None, relax=1.0, ws=None, scheme='CDS', **kwargs)[source]

Solve partial differential equation with Finite Volumes.

This function is a syntactic sugar proxy for using the Finite Volume functionality of the library core to solve elliptic and parabolic partial differential of the following type:

\[\begin{split}\frac{\partial u}{\partial t} + \mathbf{v}\cdot\nabla u & = \nabla\cdot(a \nabla u) + b u + f(\mathbf{r},t)\\ u(\mathbf{r}, t) & = u_B \quad\mathbf{r}\in\Gamma_{\text{Dirichlet}}\\ \frac{\partial u(\mathbf{r}, t)}{\partial \mathbf{n}} & = u_{\partial \text{B}} \quad\mathbf{r}\in\Gamma_{\text{Neumann}}\\ u(\mathbf{r}, t=0) & = u_0 \quad\text{with} \quad\mathbf{r}\in\Omega\end{split}\]

The Domain \(\Omega\) and the Boundary \(\Gamma\) are defined through the given mesh with appropriate boundary marker.

The solution \(u(\mathbf{r}, t)\) is given for each cell in the mesh.

Parameters
  • mesh (GIMLI::Mesh) – Mesh represents spatial discretization of the calculation domain

  • a (value | array | callable(cell, userData)) – Stiffness weighting per cell values.

  • b (value | array | callable(cell, userData)) – Scale for mass values b

  • f (iterable(cell)) – Load vector

  • fn (iterable(cell)) – TODO What is fn

  • vel (ndarray (N,dim) | RMatrix(N,dim)) – Velocity field \(\mathbf{v}(\mathbf{r}, t=\text{const}) = \{[v_i]_j,\}\) with \(i=[1\ldots 3]\) for the mesh dimension and \(j = [0\ldots N-1]\) with N either the amount of cells, nodes, or boundaries. Velocities per boundary are preferred and will be interpolated on demand.

  • u0 (value | array | callable(cell, userData)) – Starting field

  • times (iterable) – Time steps to calculate for.

  • Workspace (ws) –

    This can be an empty class that will used as an Workspace to store and cache data.

    If ws is given: The system matrix is taken from ws or calculated once and stored in ws for further usage.

    The system matrix is cached in this Workspace as ws.S The LinearSolver with the factorized matrix is cached in this Workspace as ws.solver The rhs vector is only stored in this Workspace as ws.rhs

  • scheme (str [CDS]) – Finite volume scheme: pygimli.solver.diffusionConvectionKernel

  • **kwargs

    • bc : Boundary Conditions dictionary, see pg.solver

    • uB : Dirichlet boundary conditions DEPRECATED

    • duB : Neumann boundary conditions DEPRECATED

Returns

u – Solution field for all time steps.

Return type

ndarray(nTimes, nCells)

Examples using pygimli.solver.solveFiniteVolume

triDiagToeplitz

pygimli.solver.triDiagToeplitz(dom, a, l, r, start=0, end=- 1)[source]

Create tri-diagonal Toeplitz matrix.

Classes

RungeKutta

class pygimli.solver.RungeKutta(solver, verbose=False)[source]

Bases: object

TODO DOCUMENT ME

__init__(solver, verbose=False)[source]

TODO DOCUMENT_ME

rk4a = [0.0, -0.41789047449985195, -1.192151694642677, -1.6977846924715279, -1.5141834442571558]
rk4b = [0.14965902199922912, 0.37921031299962726, 0.8229550293869817, 0.6994504559491221, 0.15305724796815198]
rk4c = [0.0, 0.14965902199922912, 0.37040095736420475, 0.6222557631344432, 0.9582821306746903]
run(u0, dt, tMax=1)[source]

TODO DOCUMENT_ME

start(u0, dt, tMax=1)[source]

TODO DOCUMENT_ME

step()[source]

TODO DOCUMENT ME

WorkSpace

class pygimli.solver.WorkSpace[source]

Bases: object