Hello,
I am trying to solve a set of coupled PDEs using MethodOfLines.jl, but I am new to both PDE solving and Julia in general. The code is based on one of the MethodOfLines.jl tutorials:
using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, DomainSets
N = 2
@parameters x t
@variables ϕ(..)[1:N]
Dt = Differential(t)
Dx = Differential(x)
Dtt = Differential(t)^2
Dxx = Differential(x)^2
γ = 0.01
α = 0.05
S = -0.499997
l = 5
bc = 0
x_min = t_min = 0
x_max = l
t_max = 20
J(ϕ) = Dtt(ϕ) + α*Dt(ϕ) - γ
eq1 = [Dxx(ϕ(x,t)[1]) ~ J(ϕ(x,t)[1]) + S*J(ϕ(x,t)[2])]
eqmid = [Dxx(ϕ(x,t)[i]) ~ S*J(ϕ(x,t)[i-1]) + J(ϕ(x,t)[i]) + S*J(ϕ(x,t)[i+1]) for i in 2:(N-1)]
eqN = [Dxx(ϕ(x,t)[N]) ~ S*J(ϕ(x,t)[N-1]) + J(ϕ(x,t)[N])]
eq = [eq1, eqmid, eqN]
eq = reduce(vcat,eq)
domains = [x ∈ Interval(x_min, x_max),
t ∈ Interval(t_min, t_max)]
bcs = [[Dx(ϕ(x_min,t)[i]) ~ bc,
Dx(ϕ(x_max,t)[i]) ~ bc,
ϕ(x,t_min)[i] ~ bc,
Dt(ϕ(x,t_min)[i]) ~ bc] for i in 1:N]
bcs = reduce(vcat,bcs)
@named pdesys = PDESystem(eq,bcs,domains,[x,t],[ϕ(x,t)[i] for i in 1:N])
Nd = 8
order = 2
discretization = MOLFiniteDifference([x=>Nd],t,approx_order=order)
println("Discretization:")
@time prob = discretize(pdesys,discretization)
This isn’t the full problem, but it results in the same error as the full problem. The error report is as follows:
ArgumentError: SymbolicUtils.BasicSymbolic{Real}[var"var\"ϕ_Any[1]\"(t)[2]ˍtt", var"var\"ϕ_Any[1]\"(t)[3]ˍtt", var"var\"ϕ_Any[1]\"(t)[4]ˍtt", var"var\"ϕ_Any[1]\"(t)[5]ˍtt", var"var\"ϕ_Any[1]\"(t)[6]ˍtt", var"var\"ϕ_Any[1]\"(t)[7]ˍtt"] are missing from the variable map.
Stacktrace:
[1] throw_missingvars(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/variables.jl:122
[2] _varmap_to_vars(varmap::Dict{Any, Any}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}}; defaults::Dict{Any, Any}, check::Bool, toterm::typeof(ModelingToolkit.default_toterm))
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/variables.jl:116
[3] varmap_to_vars(varmap::Vector{Pair}, varlist::Vector{SymbolicUtils.BasicSymbolic{Real}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/variables.jl:85
[4] get_u0_p(sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; use_union::Bool, tofloat::Bool, symbolic_u0::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:787
[5] get_u0_p
@ ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:769 [inlined]
[6] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), kwargs::@Kwargs{t::Int64, has_difference::Bool, check_length::Bool})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:812
[7] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:937
[8] ODEProblem
@ ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:930 [inlined]
[9] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:930
[10] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:917
[11] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/abstractodesystem.jl:916
...
@ PDEBase ~/.julia/packages/PDEBase/Aqj4G/src/discretization_state.jl:55
[16] macro expansion
@ ./timing.jl:279 [inlined]
Based on other posts I’ve read, the issue has to do with the boundary conditions (or lack thereof) but I can’t figure out how to get the problem to discretize even with additional BCs (the conditions at t=0 aren’t part of the problem statement). It seems to work fine if I simply don’t couple the two equations, but that honestly confuses me more.
Any help would be greatly appreciated!