Boundary Conditions

What follows is a set of allowable boundary conditions, please note that this is not exhaustive - try your condition and see if it works, the handling is quite general. If it doesn't, please post an issue and we'll try to support it. Currently, boundary conditions have to be supplied at the edge of the domain, but there are plans to support conditions embedded in the domain.

Definitions

using ModelingToolkit, MethodOfLines, DomainSets

@parameters x y t
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

x_min = y_min = 0.0

x_max = y_max = 1.0

Dirichlet

v(t, 0, y) ~ 1.0

Time dependent

u(t, 0.0, y) ~ x_min * y + 0.5t

Julia function

v(t, x, y_max) ~ sin(x)

User-defined function

alpha = 9

f(t, x, y) = x * y - t

function g(x, y)
    z = sin(x * y) + cos(y)
    # Note that symbolic conditionals require the use of Base.ifelse, or registration
    return ifelse(z > 0, x, 1.0)
end

u(t, x, y_min) ~ f(t, x, y_min) + alpha / g(x, y_min)

Registered User Defined Function

alpha = 9

f(t, x, y) = x * y - t

function g(x, y)
    z = sin(x * y) + cos(y)
    # This function must be registered as it contains a symbolic conditional
    if z > 0
        return x
    else
        return 1.0
    end
end

@register_symbolic g(x, y)

u(t, x, y_min) ~ f(t, x, y_min) + alpha / g(x, y_min)

Neumann/Robin

v(t, x_min, y) ~ 2.0 * Dx(v(t, x_min, y))

Time dependent

u(t, x_min, y) ~ x_min * Dy(v(t, x_min, y)) + 0.5t

Higher order

v(t, x, 1.0) ~ sin(x) + Dyy(v(t, x, y_max))

Time derivative

Dt(u(t, x_min, y)) ~ 0.2

User-defined function

function f(u, v)
    (u + Dyy(v) - Dy(u)) / (1 + v)
end

Dyy(u(t, x, y_min)) ~ f(u(t, x, y_min), v(t, x, y_min)) + 1

0 lhs

0 ~ u(t, x, y_max) - Dy(v(t, x, y_max))

Periodic

u(t, x_min, y) ~ u(t, x_max, y)

v(t, x, y_max) ~ u(t, x_max, y)

Please note that if you want to use a periodic condition on a dimension with WENO schemes, please use a periodic condition on all variables in that dimension.

Periodic conditions are also supported with WENOScheme() on nonuniform grids; wrapped stencils use the exact periodically shifted grid coordinates.

Interfaces

You may want to connect regions with differing dynamics together, to do this follow the following example, splitting the variable that spans these domains:

@parameters t x1 x2
@variables c1(..)
@variables c2(..)
Dt = Differential(t)

Dx1 = Differential(x1)
Dxx1 = Dx1^2

Dx2 = Differential(x2)
Dxx2 = Dx2^2

D1(c) = 1 + c / 10
D2(c) = 1 / 10 + c / 10

eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))),
    Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))]

bcs = [c1(0, x1) ~ 1 + cospi(2 * x1),
    c2(0, x2) ~ 1 + cospi(2 * x2),
    Dx1(c1(t, 0)) ~ 0,
    c1(t, 0.5) ~ c2(t, 0.5), # Relevant interface boundary condition
    -D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), # Higher order interface condition
    Dx2(c2(t, 1)) ~ 0]

domains = [t ∈ Interval(0.0, 0.15),
    x1 ∈ Interval(0.0, 0.5),
    x2 ∈ Interval(0.5, 1.0)]

@named pdesys = PDESystem(eqs, bcs, domains,
    [t, x1, x2], [c1(t, x1), c2(t, x2)])

Note that if you want to use a higher order interface condition, this may not work if you have no simple condition of the form c1(t, 0.5) ~ c2(t, 0.5).

Interfaces require the connected variables to use the same step size, except when WENOScheme() is used with nonuniform grid vectors supplied for both connected variables, in which case the grids may differ: advection stencils that cross the interface are evaluated with the exact physical coordinates of the neighbouring grid. This exception covers first-order (advection) derivatives only — systems containing higher-order spatial derivatives (e.g. diffusion) in the connected variables are rejected with an error at discretization time.