Initial and Boundary conditions

Oftentimes we will want to do a 1, 2, or 3-dimensional simulation, rather than the 0-dimensional simulation we get by default with a system of ordinary differential equations. In these cases, we will need to specify initial and boundary conditions for the system.

To demonstrate how to do this, we will use the following simple system of ordinary differential equations:

using EarthSciMLBase
using ModelingToolkit
using ModelingToolkit: t_nounits, D_nounits
t = t_nounits
D = D_nounits

@parameters x y

function ExampleSys()
    @variables u(t) v(t)
    @parameters x=1 y=2
    eqs = [
        D(u) ~ √abs(v),
        D(v) ~ √abs(u),
    ]
    ODESystem(eqs, t, [u, v], [x, y]; name=:Docs₊Example)
end

ExampleSys()

\[ \begin{align} \frac{\mathrm{d} u\left( t \right)}{\mathrm{d}t} &= \sqrt{\left|v\left( t \right)\right|} \\ \frac{\mathrm{d} v\left( t \right)}{\mathrm{d}t} &= \sqrt{\left|u\left( t \right)\right|} \end{align} \]

Next, we specify our initial and boundary conditions using the DomainInfo type. We initialize DomainInfo with sets of initial and boundary conditions. In the example below, we set constant initial conditions using constIC and constant boundary conditions using constBC.

using DomainSets

x_min = y_min = t_min = 0.0
x_max = y_max = t_max = 1.0

# Create constant initial conditions = 16.0 and boundary conditions = 4.0.
icbc = DomainInfo(
    constIC(4.0, t ∈ Interval(t_min, t_max)),
    constBC(16.0,
        x ∈ Interval(x_min, x_max),
        y ∈ Interval(y_min, y_max),
    ),
)

It is also possible to use periodic boundary conditions with periodicBC and zero-gradient boundary conditions with zerogradBC.

Finally, we combine our initial and boundary conditions with our system of equations using the couple function.

model = couple(ExampleSys(), icbc)

eq_sys = convert(PDESystem, model)

\[ \begin{align} \mathtt{{\delta}Docs.Example.x\_transform}\left( t, x, y \right) &= 1 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( t, x, y \right) &= 1 \\ \frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Docs.Example.u}\left( t, x, y \right) &= \sqrt{\left|\mathtt{Docs.Example.v}\left( t, x, y \right)\right|} \\ \frac{\mathrm{d}}{\mathrm{d}t} \mathtt{Docs.Example.v}\left( t, x, y \right) &= \sqrt{\left|\mathtt{Docs.Example.u}\left( t, x, y \right)\right|} \end{align} \]

We can also look at the expanded boundary conditions of the resulting equation system:

eq_sys.bcs

\[ \begin{align} \mathtt{{\delta}Docs.Example.x\_transform}\left( 0, x, y \right) &= 4 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( 0, x, y \right) &= 4 \\ \mathtt{Docs.Example.u}\left( 0, x, y \right) &= 4 \\ \mathtt{Docs.Example.v}\left( 0, x, y \right) &= 4 \\ \mathtt{{\delta}Docs.Example.x\_transform}\left( t, 0, y \right) &= 16 \\ \mathtt{{\delta}Docs.Example.x\_transform}\left( t, 1, y \right) &= 16 \\ \mathtt{{\delta}Docs.Example.x\_transform}\left( t, x, 0 \right) &= 16 \\ \mathtt{{\delta}Docs.Example.x\_transform}\left( t, x, 1 \right) &= 16 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( t, 0, y \right) &= 16 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( t, 1, y \right) &= 16 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( t, x, 0 \right) &= 16 \\ \mathtt{{\delta}Docs.Example.y\_transform}\left( t, x, 1 \right) &= 16 \\ \mathtt{Docs.Example.u}\left( t, 0, y \right) &= 16 \\ \mathtt{Docs.Example.u}\left( t, 1, y \right) &= 16 \\ \mathtt{Docs.Example.u}\left( t, x, 0 \right) &= 16 \\ \mathtt{Docs.Example.u}\left( t, x, 1 \right) &= 16 \\ \mathtt{Docs.Example.v}\left( t, 0, y \right) &= 16 \\ \mathtt{Docs.Example.v}\left( t, 1, y \right) &= 16 \\ \mathtt{Docs.Example.v}\left( t, x, 0 \right) &= 16 \\ \mathtt{Docs.Example.v}\left( t, x, 1 \right) &= 16 \end{align} \]