sofia package¶

Subpackages¶

Module contents¶

sofia.requires(*core_deps: str)¶: Build a decorator raising ImportError if given core dependencies are missing.

class sofia.Domain(transformation: CoordSys = global)¶

Bases: Basic

The position space where a field lives.

dim(coordinate: Symbol) → int¶: Return the dimension index of a coordinate.

vec(*components: Field)¶

property identity¶: Return a sum over the outer products of the unit vectors with themselves.

default_assumptions = {}¶

class sofia.Field(name, domain, *args, **kwargs)¶

Bases: Expr

Symbolic field living on Domain.

is_commutative = True¶

copy(new_name: str, append_name: bool = False)¶

as_dummy(name=None)¶

Return the expression with any objects having structurally bound symbols replaced with unique, canonical symbols within the object in which they appear and having only the default assumption for commutativity being True. When applied to a symbol a new symbol having only the same commutativity will be returned.

Examples

>>> from sympy import Integral, Symbol
>>> from sympy.abc import x
>>> r = Symbol('r', real=True)
>>> Integral(r, (r, x)).as_dummy()
Integral(_0, (_0, x))
>>> _.variables[0].is_real is None
True
>>> r.as_dummy()
_r

Notes

Any object that has structurally bound variables should have a property, bound_symbols that returns those symbols appearing in the object.

default_assumptions = {'commutative': True}¶

class sofia.VectorField(name: str | Symbol, domain: Domain, staggered: str | tuple[Staggered, ...] | None = None, **kwargs)¶

class sofia.VectorField(*terms: Vector, **options)

Bases: VectorAdd

A vector field composed of one scalar Field per axis.

It is a sympy vector so it can be used in any vector expression.

Parameters:

name – Each component field is named f”{name}_{axis}” (where axis is “x”, “y” or “z”).
domain – Domain the component fields live on.
staggered – Where the components sit within a cell on a Grid. This can be given either as a Staggered per component or as the string “faces”, which means each component lives on the face normal to its axis.
**kwargs – Forwarded to every component Field.

Examples

>>> from sofia import Grid
>>> G = Grid(ndim=2, name="G")
>>> u = VectorField("u", G)
>>> u
u_x*global.ihat + u_y*global.jhat
>>> u.fields
(u_x, u_y)
>>> [f.staggered for f in VectorField("m", G, staggered="faces").fields]
[XFACE, YFACE]
>>> u << 0
EvalAssignment((u_x, u_y), (0, 0))

is_commutative = True¶

property fields: tuple[Field, ...]¶: The scalar field of each component.

assignments(other) → dict[Field, Expr]¶: Return the value each component field gets when other is assigned.

default_assumptions = {'commutative': True}¶

class sofia.TensorField(name: str | Symbol, domain: Domain, symmetric: bool = False, staggered: str | tuple[tuple[Staggered | None, ...], ...] | None = None, **kwargs)¶

class sofia.TensorField(*terms: Dyadic, **options)

Bases: DyadicAdd

An order 2 tensor field (dyadic) composed of one scalar Field per pair of axes.

It is a sympy dyadic, so it can be used in any dyadic expression.

Parameters:

name – Each component field is named f”{name}_{a}{b}” (e.g. T_xy).
domain – Domain the component fields live on.
symmetric – Use one field for T_ab and T_ba, so the tensor is symmetric by construction.
staggered – Where the components sit within a cell on a Grid. This can be given either as one Staggered per component, as rows (staggered[a][b] for component T_ab), or as the string “mac”, which puts the diagonal components at cell centers and the off-diagonal ones on the edges.
**kwargs – Forwarded to every component Field.

Examples

>>> from sofia import Grid
>>> G = Grid(ndim=2, name="G")
>>> T = TensorField("T", G)
>>> T.fields
(T_xx, T_xy, T_yx, T_yy)
>>> T.dot(G.unit_vectors[0])
T_xx*global.ihat + T_yx*global.jhat
>>> S = TensorField("S", G, symmetric=True)
>>> S.fields
(S_xx, S_xy, S_yy)
>>> S << 2
EvalAssignment((S_xx, S_xy, S_yy), (2, 2, 2))

is_commutative = True¶

property fields: tuple[Field, ...]¶: The distinct component fields.

property symmetric: bool¶: Whether T_ab and T_ba are the same field for every pair of axes.

assignments(other) → dict[Field, Expr]¶: Return the value each component field gets when other is assigned.

default_assumptions = {'commutative': True}¶

Bases: Field

Symbolic representation of a field living on a structured grid.

A GridField is a symbol that carries the information needed to discretise expressions it appears in: the grid it is defined on, which of the grid axes it varies along, where within a cell its values sit, and how far it is shifted from its natural location.

Parameters:

name – Name of the field, used for its symbolic and generated-code names.
domain – The Grid the field is defined on.
staggered – Location of the field values within a cell, as a Staggered member. Defaults to Staggered.CENTER.
discretisation_order – Order of accuracy for derivatives and interpolation of this field. Must be even. The default of 0 means the order of domain is used instead.
axes – Grid axes the field varies along, given as indices or coordinates. Defaults to all of them. Passing a subset gives a field of lower dimensionality than the grid; passing () gives a scalar.

Examples

>>> import sympy as sp
>>> from sofia import Grid, Field, Staggered

>>> G = Grid(ndim=2)
>>> u = Field("u", domain=G)

>>> u.shape
(G.Nx, G.Ny)
>>> u.axes
(0, 1)
>>> Field("v", domain=G, axes=(G.x,)).shape
(G.Nx,)
>>> Field("dt", domain=G, axes=()).shape
()

Shifts:

>>> u.s(1, 0).shifts
(1, 0)
>>> u.s(1, 0).s(2, 0).shifts
(3, 0)
>>> u.s(j=-1).shifts
(0, -1)

is_commutative = True¶

copy(new_name: str, append_name: bool = False)¶

with_order(order)¶

to_staggered_origin()¶

reset_shifts()¶

s(*shifts, **kwargs)¶

indexify(bare=False)¶

fractional_offset(wrt) → Expr¶: Combine shift and staggered offset of a field.

interpolate()¶

default_assumptions = {'commutative': True}¶

Bases: Domain

Symbolic representation of a structured grid.

This class defines a discretized grid in arbitrary dimensional space using symbolic quantities. A Grid instance serves multiple purposes when building a simulation code, but most directly it can be used to conveniently get symbolic coordinates and indices, as well as symbolic grid shape and size specifications.

Parameters:

ndim – Number of dimensions of the grid.
name – Name of the grid. This name is used to generate symbolic names for the grid coordinates, indices, and other symbolic quantities.
discretisation_order – Order of the discretisation scheme. Must be an even positive integer.
boundaries – Boundary specifications for each dimension in index space. Each dimension can be specified as either a single value (for periodic boundaries) or a pair of values (for non-periodic boundaries). The values represent the insets from first and last grid nodes. For example, a boundary specification of (0.5, -2.5) means that the boundary is at i = 0.5 and i = N - 1 - 2.5, where N is the number of grid points in that dimension.
transformations – Map from computational to physical coordinates, as a tuple of callables, one per axis. It may contain a transformation placeholder, an undefined sympy function applied to the axis coordinate, which makes the axis discrete. The physical coordinate of that axis is then the field G.transformation_fields[c], with c the axis coordinate.
transformations_inverse – Map from physical back to computational coordinates, in the same form. If None, it is solved for with sp.solve. Ignored when any axis has a transformation placeholder; such a grid has no inverse.

Examples

A discrete x axis, with the transformation placeholder as the map or inside one:

>>> X = sp.Function("X")
>>> G1 = Grid(ndim=2, name="G1", transformations=(X, lambda y: y))
>>> G2 = Grid(ndim=2, name="G2", transformations=(lambda x: 2 * X(x), lambda y: y))
>>> G1.scale_factors
(Derivative(X, G1.x), 1)
>>> G2.scale_factors
(2*Derivative(X, G2.x), 1)

The physical coordinate of the discrete axis is a field on the grid, passed to the simulation as an input:

>>> X1 = G1.transformation_fields[G1.x]
>>> type(X1).__name__, X1.domain == G1
('GridField', True)

property uniform_axes¶: Axes whose scale factor is constant – these stay circulant/spectral.

property transformation_fields¶: Map each discrete axis’s coordinate to the field of its node positions.

property boundaries_absolute¶: Physical domain edges in index space.

property exterior_nodes¶: Number of nodes per dim that lie outside the physcal domain (ghost cells).

dim(coordinate_or_index: Symbol) → int¶: Return the dimension index of a coordinate.

coord_to_ind(coordinate: Symbol) → Symbol¶

dds_eval(L, boundary_left, boundary_right)¶

maps_coord_to_ind()¶

maps_ind_to_coord()¶

Return the substitutions corresponding to the concrete numbers of the grid.

Parameters:

shape – Grid shape per dimension.
extent – Physical domain extents per dimension, given as concrete floats.
subs – Any additional substitutions.

Examples

>>> from sofia import Grid

>>> G = Grid(ndim=2)
>>> G.concrete(shape=(100, 200), extent=((0, 1), (0, 2)))
{G.Nx: 100, G.Ny: 200, G.Lx: 1, G.Ly: 2, G.x0: 0, G.y0: 0, global.x: i/100, global.y: j/100}

default_assumptions = {}¶

class sofia.NeighbourSum(*args)¶

Bases: Function

Sum a per-pair expression over the neighbours of each particle.

Parameters:

expr – Per-pair expression. field is the value at the particle itself, field.neighbour the value at the neighbour. Can be scalar or vector.
pos – Particle position, determines which particles count as neighbours.
h – Cutoff radius, only pairs closer than h contribute to the sum.
exclude_self – If True, skip the particle itself. Only for sums over the same domain.
neighbours – Positions of the particles to sum over, on this or another particle domain. Defaults to pos. In a sum over another domain, fields of that domain are the values at the neighbour. .neighbour is only needed when both are on the same domain.

property cutoff: Basic¶

Condition for a neighbour to contribute to the sum.

A neighbour contributes if it is closer than h to the particle.

default_assumptions = {}¶

class sofia.ParticleDomain(ndim: Literal[1, 2, 3] = 3, name: str = 'P', _n_particles: Expr | int | None = None, _extent: Tuple | None = None)¶

Bases: Domain

Symbolic representation of a particle domain.

A domain where fields are indexed by particle. Particles live in ndim-dimensional Cartesian space.

Parameters:

ndim – Spatial dimensionality.
name – Name used to identify the particle domain.

concrete(n_particles: int | None = None, extent: tuple[Numeric | tuple[Numeric, Numeric] | list[Numeric], ...] | list[Numeric | tuple[Numeric, Numeric] | list[Numeric]] | Numeric | None = None, subs: dict | None = None) → dict[Symbol, int | float]¶: Return substitutions corresponding to the concrete numbers of the domain.

default_assumptions = {}¶

class sofia.ParticleField(name: str, domain: ParticleDomain)¶

Bases: Field

Symbolic field living on a ParticleDomain.

indexify()¶

default_assumptions = {'commutative': True}¶

is_commutative = True¶

class sofia.Scalar(name: str | Symbol, *, dtype: DTypeLike | None = None, **assumptions: bool | None)¶

Bases: AtomicExpr

Symbolic scalar.

Identified by name and assumptions, compared after sympy has drawn their consequences. If dtype is boolean, a BooleanScalar is returned instead.

Parameters:

name – Identifier the generated code declares this scalar under.
dtype – Data type. By default it is the backend’s float dtype.
assumptions – Sympy assumptions, restricted to SUPPORTED_ASSUMPTIONS.

is_commutative = True¶

is_symbol = True¶

is_Symbol = True¶

kind = NumberKind¶

default_assumptions = {'commutative': True}¶

class sofia.Staggered(offsets: dict[int, Expr | float] | None = None, name=None)¶

Bases: Basic

Fractional per-axis offsets of a field relative to cell centers.

CENTER: ClassVar[Staggered] = Center¶

XFACE: ClassVar[Staggered] = XFACE¶

YFACE: ClassVar[Staggered] = YFACE¶

ZFACE: ClassVar[Staggered] = ZFACE¶

staggered_offsets(ndim: int) → tuple[Expr, ...]¶

XY_EDGE = XY_EDGE¶

XZ_EDGE = XZ_EDGE¶

YZ_EDGE = YZ_EDGE¶

default_assumptions = {}¶

sofia.fields(names: str, domain: Domain, **kwargs: object) → tuple[Field, ...]¶

Create one or more Fields on a shared domain.

Equivalent to sympy.symbols(), but for Fields. Always returns a tuple for better type support.

Parameters:

names – Field names, comma- or whitespace-separated; see sympy.symbols().
domain – Domain the fields live on.
**kwargs – Forwarded to the Field subclass registered for domain.

Returns:

A tuple of the created Field(s).

Return type:

tuple[Field, …]

Examples

>>> (u,) = fields("u", domain=Grid(ndim=1))
>>> v, w = fields("v w", domain=Grid(ndim=2))

sofia.print_environment_info() → None¶: Print the environment info for reporting bugs, mimicking the jax equivalent.

sofia.scalars(names: str, **kwargs: object) → tuple[Scalar, ...]¶

Create one or more Scalars.

Equivalent to sympy.symbols(), but for Scalars. Always returns a tuple for better type support.

Parameters:

names – Scalar names, comma- or whitespace-separated; see sympy.symbols().
**kwargs – Keyword arguments passed to Scalar.

Returns:

A tuple of the created Scalar(s).

Return type:

tuple[Scalar, …]

Examples

>>> (a,) = scalars("a")
>>> b, c = scalars("b c", positive=True)