Allow Interface boundary conditions

docs/src/boundary_conditions.md

@@ -100,4 +100,40 @@ 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.
+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.
+
+## 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:
+```julia
+    @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)`.

src/MOL_discretization.jl

@@ -30,11 +30,14 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
     bcorders = Dict(map(x -> x => d_orders(x, pdesys.bcs), all_ivs(v)))
     # Create a map of each variable to their boundary conditions including initial conditions
     boundarymap = parse_bcs(pdesys.bcs, v, bcorders)
-    # Generate a map of each variable to whether it is periodic in a given direction
-    pmap = PeriodicMap(boundarymap, v)
+
     # Transform system so that it is compatible with the discretization
     if discretization.should_transform
-        pdesys, pmap = transform_pde_system!(v, boundarymap, pmap, pdesys)
+        if has_interfaces(boundarymap)
+            @warn "The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature."
+        else
+            pdesys = transform_pde_system!(v, boundarymap, pdesys)
+        end
     end
 
     # Check if the boundaries warrant using ODAEProblem, as long as this is allowed in the interface
@@ -67,7 +70,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
     s = DiscreteSpace(v, discretization)
     # Get the interior and variable to solve for each equation
     #TODO: do the interiormap before and independent of the discretization i.e. `s`
-    interiormap = InteriorMap(pdeeqs, boundarymap, s, discretization, pmap)
+    interiormap = InteriorMap(pdeeqs, boundarymap, s, discretization)
     # Get the derivative orders appearing in each equation
     pdeorders = Dict(map(x -> x => d_orders(x, pdeeqs), v.x̄))
     bcorders = Dict(map(x -> x => d_orders(x, bcs), v.x̄))
@@ -117,7 +120,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
             if disc_strategy isa ScalarizedDiscretization
                 # Generate the equations for the interior points
                 discretize_equation!(alleqs, bceqs, pde, interiormap, eqvar, bcmap,
-                    depvars, s, derivweights, indexmap, pmap)
+                    depvars, s, derivweights, indexmap)
             else
                 throw(ArgumentError("Only ScalarizedDiscretization is currently supported"))
             end
@@ -139,15 +142,15 @@ function SciMLBase.discretize(pdesys::PDESystem,discretization::MethodOfLines.MO
     try
         simpsys = structural_simplify(sys)
         if tspan === nothing
-            add_metadata!(simpsys.metadata, sys)
+            add_metadata!(get_metadata(sys), sys)
             return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); discretization.kwargs...)
         else
             # Use ODAE if nessesary
             if getfield(sys, :metadata) isa MOLMetadata && getfield(sys, :metadata).use_ODAE
-                add_metadata!(simpsys.metadata, DAEProblem(simpsys; discretization.kwargs...))
+                add_metadata!(get_metadata(simpsys), DAEProblem(simpsys; discretization.kwargs...))
                 return prob = ODAEProblem(simpsys, Pair[], tspan; discretization.kwargs...)
             else
-                add_metadata!(simpsys.metadata, sys)
+                add_metadata!(get_metadata(simpsys), sys)
                 return prob = ODEProblem(simpsys, Pair[], tspan; discretization.kwargs...)
             end
         end
@@ -158,6 +161,7 @@ end
 
 function get_discrete(pdesys, discretization)
     t = discretization.time
+    disc_strategy = discretization.disc_strategy
     cardinalize_eqs!(pdesys)
 
     ############################
@@ -173,11 +177,14 @@ function get_discrete(pdesys, discretization)
     bcorders = Dict(map(x -> x => d_orders(x, pdesys.bcs), all_ivs(v)))
     # Create a map of each variable to their boundary conditions including initial conditions
     boundarymap = parse_bcs(pdesys.bcs, v, bcorders)
-    # Generate a map of each variable to whether it is periodic in a given direction
-    pmap = PeriodicMap(boundarymap, v)
+
     # Transform system so that it is compatible with the discretization
     if discretization.should_transform
-        pdesys, pmap = transform_pde_system!(v, boundarymap, pmap, pdesys)
+        if has_interfaces(boundarymap)
+            @warn "The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature."
+        else
+            pdesys = transform_pde_system!(v, boundarymap, pdesys)
+        end
     end
 
     s = DiscreteSpace(v, discretization)

src/MOL_utils.jl

@@ -269,29 +269,6 @@ remove(v::AbstractVector, a::Number) = filter(x -> !isequal(x, a), v)
 
 half_range(x) = -div(x, 2):div(x, 2)
 
-@inline function _wrapperiodic(I, N, j, l)
-    I1 = unitindex(N, j)
-    # -1 because of the relation u[1] ~ u[end]
-    if I[j] <= 1
-        I = I + I1 * (l - 1)
-    elseif I[j] > l
-        I = I - I1 * (l - 1)
-    end
-    return I
-end
-
-"""
-Allow stencils indexing over periodic boundaries. Index through this function.
-"""
-@inline function wrapperiodic(I, s, ::Val{true}, u, jx)
-    j, x = jx
-    return _wrapperiodic(I, ndims(u, s), j, length(s, x))
-end
-
-@inline function wrapperiodic(I, s, ::Val{false}, u, jx)
-    return I
-end
-
 d_orders(x, pdeeqs) = reverse(sort(collect(union((differential_order(pde.rhs, safe_unwrap(x)) for pde in pdeeqs)..., (differential_order(pde.lhs, safe_unwrap(x)) for pde in pdeeqs)...))))
 
 insert(args...) = insert!(args[1], args[2:end]...)

src/MethodOfLines.jl

@@ -3,7 +3,7 @@ using LinearAlgebra
 using SciMLBase
 using DiffEqBase
 using ModelingToolkit
-using ModelingToolkit: operation, istree, arguments, variable
+using ModelingToolkit: operation, istree, arguments, variable, get_metadata
 using SymbolicUtils, Symbolics
 using Symbolics: unwrap, solve_for, expand_derivatives, diff2term, setname, rename, similarterm
 using SymbolicUtils: operation, arguments
@@ -16,7 +16,12 @@ import DomainSets
 # To Extend
 import SciMLBase.wrap_sol
 import Base.display
-
+import Base.isequal
+import Base.getindex
+import Base.checkindex
+import Base.checkbounds
+import Base.getproperty
+import Base.ndims
 
 # Interface
 include("interface/grid_types.jl")
@@ -50,6 +55,9 @@ include("system_parsing/bcs/parse_boundaries.jl")
 include("system_parsing/bcs/periodic_map.jl")
 include("system_parsing/pde_system_transformation.jl")
 
+# Interface handling
+include("discretization/interface_boundary.jl")
+
 # Schemes
 include("discretization/schemes/centered_difference/centered_difference.jl")
 include("discretization/schemes/upwind_difference/upwind_difference.jl")

src/discretization/discretize_vars.jl

@@ -101,7 +101,6 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
         end
         x => discx
     end
-
     # Define the grid on which the dependent variables will be evaluated (see #378)
     # center_align is recommended for Dirichlet BCs
     # edge_align is recommended for Neumann BCs (spatial discretization is conservative)
@@ -149,7 +148,6 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
     return DiscreteSpace{nspace,length(depvars),G}(vars, Dict(depvarsdisc), axies, grid, Dict(dxs), Dict(Iaxies), Dict(Igrid))
 end
 
-import Base.getproperty
 
 function Base.getproperty(s::DiscreteSpace, p::Symbol)
     if p in [:ū, :x̄, :time, :args, :x2i, :i2x]
@@ -172,7 +170,7 @@ nvars(::DiscreteSpace{N,M}) where {N,M} = M
 Fillter out the time variable and get the spatial variables of `u` in `s`.
 """
 params(u, s::DiscreteSpace) = remove(s.args[operation(u)], s.time)
-Base.ndims(u, s::DiscreteSpace) = length(params(u, s))
+Base.ndims(u, s::DiscreteSpace) = ndims(s.discvars[depvar(u, s)])
 
 Base.length(s::DiscreteSpace, x) = length(s.grid[x])
 Base.length(s::DiscreteSpace, j::Int) = length(s.grid[s.x̄[j]])
@@ -187,6 +185,7 @@ of `II` that corresponds to only the spatial arguments of `u`.
 @inline function Idx(II::CartesianIndex, s::DiscreteSpace, u, indexmap)
     # We need to construct a new index as indices may be of different size
     length(params(u, s)) == 0 && return CartesianIndex()
+    !all(x -> haskey(indexmap, x), params(u, s)) && return II
     is = [II[indexmap[x]] for x in params(u, s)]
 
     II = CartesianIndex(is...)