Now updated to MethodOfLines@0.11.3 (i.e. yesterday’s release), adapted my code according to updated manual - it works now 
Also got the Tutorial example with an interface working (just had to change the algorithm in my code to Rosenbrock23).
Now trying to adapt that code with interface to my 2D problem. The code below describes a simplified case of a heat/mass exchanger with a counterflow, and semi-permeable wall in the middle.
Code counterflow exchanger
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets
@parameters t x y1 y2
@variables c1(..)
@variables c2(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Dx^2
Dy1 = Differential(y1)
Dyy1 = Dy1^2
Dy2 = Differential(y2)
Dyy2 = Dy2^2
li_vel(y) = 2 * (0.5 - y)
eq1 = Dt(c1(t, x, y1)) ~ Dxx(c1(t, x, y1)) + Dyy1(c1(t, x, y1)) - Dx(c1(t, x, y1))*li_vel(y1)
eq2 = Dt(c2(t, x, y2)) ~ Dxx(c2(t, x, y2)) + Dyy2(c2(t, x, y2)) - Dx(c2(t, x, y2))*li_vel(y2)
eqs = [eq1, eq2]
domains = [t ∈ Interval(0.0, 0.15),
x ∈ Interval(0.0, 1.0),
y1 ∈ Interval(0.0, 0.5),
y2 ∈ Interval(0.5, 1.0)]
bcs = [
c1(0, x, y1) ~ 0.5, # init
c2(0, x, y2) ~ 0.5,
c1(t, 0, y1) ~ 0, # flow in
c2(t, 1, y2) ~ 1,
Dx(c1(t, 1, y1)) ~ 0, # flow out
Dx(c2(t, 0, y1)) ~ 0,
Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
Dy2(c2(t, x, 1)) ~ 0,
Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5),
]
domains = [t ∈ Interval(0.0, 0.15),
y1 ∈ Interval(0.0, 0.5),
y2 ∈ Interval(0.5, 1.0)]
@named pdesys = PDESystem(eqs, bcs, domains,
[t, x, y1, y2], [c1(t, x, y1), c2(t, x, y2)])
discretization = MOLFiniteDifference([x => 0.1, y1 => 0.1, y2 => 0.1], t)
prob = MethodOfLines.discretize(pdesys, discretization)
# error on previous line
alg = Rosenbrock23()
sol = solve(prob, alg, saveat = 0.1);
Unfortunately I get an error on MethodOfLines.discretize(pdesys, discretization):
ERROR: ArgumentError: invalid index: nothing of type Nothing
Stack trace
ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
[1] to_index(i::Nothing)
@ Base ./indices.jl:300
[2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing)
@ Base ./indices.jl:277
[3] _to_indices1(A::Vector{Symbolics.VarDomainPairing}, inds::Tuple{Base.OneTo{Int64}}, I1::Nothing)
@ Base ./indices.jl:359
[4] to_indices
@ ./indices.jl:354 [inlined]
[5] to_indices
@ ./indices.jl:345 [inlined]
[6] getindex
@ ./abstractarray.jl:1291 [inlined]
[7] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real})
@ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:34
[8] iterate
@ ./generator.jl:47 [inlined]
[9] collect_to!(dest::Vector{Pair{…}}, itr::Base.Generator{Vector{…}, PDEBase.var"#54#62"{…}}, offs::Int64, st::Int64)
@ Base ./array.jl:892
[10] collect_to_with_first!(dest::Vector{…}, v1::Pair{…}, itr::Base.Generator{…}, st::Int64)
@ Base ./array.jl:870
[11] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
@ Base ./array.jl:864
[12] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#54#62"{Vector{…}}})
@ Base ./array.jl:763
[13] map(f::Function, A::Vector{Any})
@ Base ./abstractarray.jl:3285
[14] PDEBase.VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference{…}; replaced_vars::Dict{…})
@ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:33
[15] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/symbolic_discretize.jl:17
[16] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:57
[17] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:54
[18] top-level scope
@ ~/Julia/MembH2Otransport.jl/src/PDE/m20_2D_middle-interf_4discourse.jl:52
Some type information was truncated. Use `show(err)` to see complete types.
Any ideas?
The manual says
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)
I do not have that kind of conditions. However if I add a simple condition,
bcs with a simple condition
bcs = [
c1(0, x, y1) ~ 0.5, # init
c2(0, x, y2) ~ 0.5,
c1(t, 0, y1) ~ 0, # flow in
c2(t, 1, y2) ~ 1,
Dx(c1(t, 1, y1)) ~ 0, # flow out
Dx(c2(t, 0, y1)) ~ 0,
Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
Dy2(c2(t, x, 1)) ~ 0,
Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
# Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5),
c2(t, x, 0.5) ~ c1(t, x, 0.5), # simple condition
]
I still get the same error, so it is probably not the error cause, but the next problem.