Quick start

SymPDE builds symbolic differential operators, variational forms, and PDE problems on explicitly declared domains and function spaces. User-facing objects are available from sympde.api.

Domains and fields

A field belongs to a scalar or vector function space. The following complete example creates two scalar fields on a square and forms a differential expression from them:

domain = Square('Omega')
space = ScalarFunctionSpace('V', domain)
trial, test = elements_of(space, names='u, v')

strong_expression = -div(grad(trial)) + trial
assert strong_expression.has(trial)

Differential operators remain symbolic, while linearity and product rules are applied automatically. SymPDE also supports curl, laplace, normal derivatives, traces, and interface operators.

Variational forms

Integrals identify the domain on which an expression is evaluated. A BilinearForm declares its trial and test fields, and a LinearForm declares its test field:

bilinear_form = BilinearForm(
    (trial, test),
    integral(domain, dot(grad(trial), grad(test)) + trial * test),
)

x, y = domain.coordinates
linear_form = LinearForm(test, integral(domain, cos(x + y) * test))

equation = find(
    trial,
    forall=test,
    lhs=bilinear_form(trial, test),
    rhs=linear_form(test),
)

The sympde.api.find() helper combines those forms into a symbolic equation. It is possible to impose essential boundary conditions (BCs) on the solution (typically, homogeneous Dirichlet BCs) by passing an EssentialBC object through the optional bc argument.

Expression lowering

Generic vector-calculus operators can be lowered to coordinate derivatives for a specific domain with TerminalExpr:

coordinate_expression = TerminalExpr(
    dot(grad(trial), grad(test)),
    domain,
)
assert coordinate_expression.has(trial, test)

Expression trees

Like any SymPy expression, a SymPDE expression is a tree of operators and operands. SymPy’s dotprint function writes this tree in the DOT language of Graphviz. For the operator \(\mathcal{L} u := - \nabla^2 u + u\):

def L(w):
    return -div(grad(w)) + w


operator_graph = dotprint(L(u))

digraph{ # Graph style "ordering"="out" "rankdir"="TD" ######### # Nodes # ######### "Add(Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True)))), ScalarFunction('u', commutative=True))_()" ["color"="black", "label"="Add", "shape"="ellipse"]; "Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True))))_(0,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "Integer(-1)_(0, 0)" ["color"="black", "label"="-1", "shape"="ellipse"]; "Div(Grad(ScalarFunction('u', commutative=True)))_(0, 1)" ["color"="black", "label"="Div", "shape"="ellipse"]; "Grad(ScalarFunction('u', commutative=True))_(0, 1, 0)" ["color"="black", "label"="Grad", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(0, 1, 0, 0)" ["color"="black", "label"="u", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(1,)" ["color"="black", "label"="u", "shape"="ellipse"]; ######### # Edges # ######### "Add(Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True)))), ScalarFunction('u', commutative=True))_()" -> "Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True))))_(0,)"; "Add(Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True)))), ScalarFunction('u', commutative=True))_()" -> "ScalarFunction('u', commutative=True)_(1,)"; "Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True))))_(0,)" -> "Integer(-1)_(0, 0)"; "Mul(Integer(-1), Div(Grad(ScalarFunction('u', commutative=True))))_(0,)" -> "Div(Grad(ScalarFunction('u', commutative=True)))_(0, 1)"; "Div(Grad(ScalarFunction('u', commutative=True)))_(0, 1)" -> "Grad(ScalarFunction('u', commutative=True))_(0, 1, 0)"; "Grad(ScalarFunction('u', commutative=True))_(0, 1, 0)" -> "ScalarFunction('u', commutative=True)_(0, 1, 0, 0)"; }

Expression tree of \(\mathcal{L} u := - \nabla^2 u + u\).

Generic operators such as grad and div are kept unevaluated. The atomic derivatives dx and dy are linear operators that apply the product rule immediately, so the tree of \(\partial_x \partial_y (uv)\) contains the expanded result:

expr = dx(dy(u * v))
assert expr == u * dx(dy(v)) + v * dx(dy(u)) + dx(u) * dy(v) + dx(v) * dy(u)

derivative_graph = dotprint(expr)

digraph{ # Graph style "ordering"="out" "rankdir"="TD" ######### # Nodes # ######### "Add(Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True)))), Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True)))), Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True))), Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True))))_()" ["color"="black", "label"="Add", "shape"="ellipse"]; "Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True))))_(0,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(0, 0)" ["color"="black", "label"="u", "shape"="ellipse"]; "dx(dy(ScalarFunction('v', commutative=True)))_(0, 1)" ["color"="black", "label"="dx", "shape"="ellipse"]; "dy(ScalarFunction('v', commutative=True))_(0, 1, 0)" ["color"="black", "label"="dy", "shape"="ellipse"]; "ScalarFunction('v', commutative=True)_(0, 1, 0, 0)" ["color"="black", "label"="v", "shape"="ellipse"]; "Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True))))_(1,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "ScalarFunction('v', commutative=True)_(1, 0)" ["color"="black", "label"="v", "shape"="ellipse"]; "dx(dy(ScalarFunction('u', commutative=True)))_(1, 1)" ["color"="black", "label"="dx", "shape"="ellipse"]; "dy(ScalarFunction('u', commutative=True))_(1, 1, 0)" ["color"="black", "label"="dy", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(1, 1, 0, 0)" ["color"="black", "label"="u", "shape"="ellipse"]; "Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True)))_(2,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "dx(ScalarFunction('u', commutative=True))_(2, 0)" ["color"="black", "label"="dx", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(2, 0, 0)" ["color"="black", "label"="u", "shape"="ellipse"]; "dy(ScalarFunction('v', commutative=True))_(2, 1)" ["color"="black", "label"="dy", "shape"="ellipse"]; "ScalarFunction('v', commutative=True)_(2, 1, 0)" ["color"="black", "label"="v", "shape"="ellipse"]; "Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True)))_(3,)" ["color"="black", "label"="Mul", "shape"="ellipse"]; "dx(ScalarFunction('v', commutative=True))_(3, 0)" ["color"="black", "label"="dx", "shape"="ellipse"]; "ScalarFunction('v', commutative=True)_(3, 0, 0)" ["color"="black", "label"="v", "shape"="ellipse"]; "dy(ScalarFunction('u', commutative=True))_(3, 1)" ["color"="black", "label"="dy", "shape"="ellipse"]; "ScalarFunction('u', commutative=True)_(3, 1, 0)" ["color"="black", "label"="u", "shape"="ellipse"]; ######### # Edges # ######### "Add(Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True)))), Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True)))), Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True))), Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True))))_()" -> "Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True))))_(0,)"; "Add(Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True)))), Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True)))), Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True))), Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True))))_()" -> "Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True))))_(1,)"; "Add(Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True)))), Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True)))), Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True))), Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True))))_()" -> "Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True)))_(2,)"; "Add(Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True)))), Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True)))), Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True))), Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True))))_()" -> "Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True)))_(3,)"; "Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True))))_(0,)" -> "ScalarFunction('u', commutative=True)_(0, 0)"; "Mul(ScalarFunction('u', commutative=True), dx(dy(ScalarFunction('v', commutative=True))))_(0,)" -> "dx(dy(ScalarFunction('v', commutative=True)))_(0, 1)"; "dx(dy(ScalarFunction('v', commutative=True)))_(0, 1)" -> "dy(ScalarFunction('v', commutative=True))_(0, 1, 0)"; "dy(ScalarFunction('v', commutative=True))_(0, 1, 0)" -> "ScalarFunction('v', commutative=True)_(0, 1, 0, 0)"; "Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True))))_(1,)" -> "ScalarFunction('v', commutative=True)_(1, 0)"; "Mul(ScalarFunction('v', commutative=True), dx(dy(ScalarFunction('u', commutative=True))))_(1,)" -> "dx(dy(ScalarFunction('u', commutative=True)))_(1, 1)"; "dx(dy(ScalarFunction('u', commutative=True)))_(1, 1)" -> "dy(ScalarFunction('u', commutative=True))_(1, 1, 0)"; "dy(ScalarFunction('u', commutative=True))_(1, 1, 0)" -> "ScalarFunction('u', commutative=True)_(1, 1, 0, 0)"; "Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True)))_(2,)" -> "dx(ScalarFunction('u', commutative=True))_(2, 0)"; "Mul(dx(ScalarFunction('u', commutative=True)), dy(ScalarFunction('v', commutative=True)))_(2,)" -> "dy(ScalarFunction('v', commutative=True))_(2, 1)"; "dx(ScalarFunction('u', commutative=True))_(2, 0)" -> "ScalarFunction('u', commutative=True)_(2, 0, 0)"; "dy(ScalarFunction('v', commutative=True))_(2, 1)" -> "ScalarFunction('v', commutative=True)_(2, 1, 0)"; "Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True)))_(3,)" -> "dx(ScalarFunction('v', commutative=True))_(3, 0)"; "Mul(dx(ScalarFunction('v', commutative=True)), dy(ScalarFunction('u', commutative=True)))_(3,)" -> "dy(ScalarFunction('u', commutative=True))_(3, 1)"; "dx(ScalarFunction('v', commutative=True))_(3, 0)" -> "ScalarFunction('v', commutative=True)_(3, 0, 0)"; "dy(ScalarFunction('u', commutative=True))_(3, 1)" -> "ScalarFunction('u', commutative=True)_(3, 1, 0)"; }

Expression tree of \(\partial_x \partial_y (uv)\).

To draw a tree outside of the documentation, save the output of dotprint to a file such as tree.dot and run dot -Tpng tree.dot -o tree.png.

Mappings and multipatch domains

Mappings turn logical domains into physical domains. Multiple patches can be joined by declaring their connected sides and orientation. See the Examples for analytical mappings, interface forms, exterior calculus, and the linearization of nonlinear problems. Each example is a complete script that the documentation workflow runs before building the pages.

Printing

The public latex function produces LaTeX for expressions and forms:

latex_equation = latex(equation)
print(latex_equation)

The documentation workflow executes quickstart.py before building these pages, so its imports and demonstrated API are checked continuously.