MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

Published 20 hours ago • verified publisher icon SciML

Differential w.r.t. multiple variables Any[] not allowed during MOLFiniteDifference discretize

Open BradyPlanden opened this issue 2 years ago • 6 comments

I'm trying to discretize the following Nerst-Plank system in a single spatial dimension; however, I've received the following error during the check_equations() call:

ERROR: ArgumentError: Differential w.r.t. multiple variables Any[τ, Z] are not allowed.

Full code:


using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters Z τ
@variables c(..) ϕ(..) S(..)
Dz = Differential(Z)
Dzz = Differential(Z)^2
Dτ = Differential(τ)

Z_min = Z_min = τ_min = 0.0
Z_max = 1.0
τ_max = 50.0

# Parameters
iₐ = 10.
i₀ = 20.
D₁ = 1e-11
F = 96487.
R = 8.314
T = 298.15
Lₓ = 50e-6
c₀ = 1000.
Mᵥ = 1.2998e-5
Cᵣ = cᵥ = 1.

# Dimensionless Parameters
δ = iₐ*Lₓ*Cᵣ/(2*F*D₁*c₀)
tD = abs(Cᵣ)*Lₓ^2/(3600*D₁)
Da = i₀*Lₓ/(2*F*D₁*c₀)
Pe = Mᵥ*i₀*3600/(Lₓ*Cᵣ*F)
ϕ₀ = R*T/F

# PDE
eq = [
       Dz(c(Z,τ))*Dz(ϕ(Z,τ))+c(Z,τ)*(Dzz(ϕ(Z,τ)))/(1.0-S(τ))^2 ~ 0.0, 
       Dτ(c(Z,τ)) ~ (-Pe*(c(Z,τ)/c₀)^(1/2)*(-ϕ(Z,τ)/ϕ₀)*(1.0-Z)*(Dz(c(Z,τ))))/(tD*(1.0-S(τ)))+(Dzz(c(Z,τ))/(tD*(1.0-S(τ))^2))
]
    
    
domains = [Z ∈ Interval(Z_min, Z_max),
           τ ∈ Interval(τ_min, τ_max)]

# BCs
bcs =  [   
        c(Z,0) ~ 1.0, #Concentration @ τ=0
        c(Z_max,τ)*(Dz(ϕ(1,τ))) ~ -δ*(1-S(τ)), #ϕ @ Z = 1
        ϕ(Z_max,τ) ~ 0.0, #ϕ @ Z = 1
        Dz(c(Z_max,τ)) ~ -δ*(1-S(τ)), #Concentration @ Z = 1
        Dz(ϕ(Z_max,τ)) ~ 0.0, #ϕ @ Z = 0
        Dz(ϕ(Z_min,τ))/(1-S(τ)) ~ -Da*(c(Z_min,τ)^(1/2))-ϕ(Z_min,τ), #ϕ BC @ Z = 0
        Dz(c(Z_min,τ))/(1-S(τ)) ~ -Da*c(Z_min,τ)^(1/2)*(-ϕ(Z_min,τ))+cᵥ*Pe*c(Z_min,τ)^(3/2)*(-ϕ(Z_min,τ)) #Concentration @ Z = 0
       ] 

@named pdesys = PDESystem(eq,bcs,domains,[Z,τ],[c(Z,τ), ϕ(Z,τ)])

N = 10

dz = (Z_max-Z_min)/N

order = 2

discretization = MOLFiniteDifference([Z => dz], τ, approx_order=order, grid_align=center_align)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

Full error:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Differential w.r.t. multiple variables Any[τ, Z] are not allowed.
Stacktrace:
 [1] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:171
 [2] ODESystem(deqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}}, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Any}, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, connections::Nothing, preface::Nothing, events::Vector{ModelingToolkit.SymbolicContinuousCallback}, tearing_state::Nothing, substitutions::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:120
 [3] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:174
 [4] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [6] top-level scope
   @ ~/Documents/Git/test/Scratch1.jl:143

S(τ) is a moving boundary with velocity defined as:

Dτ(S(τ)) ~ (Pe*ϕ(Z,τ)*c(Z,τ)^(1/2))/(τ*D₁)

that I am planning to combine and solve with the nerst-plank discretisation. If anyone has a better approach to this, I'm definitely open to suggestions.

BradyPlanden avatar Apr 13 '22 16:04 BradyPlanden

You'll want to include your equation for Dτ(S(τ)) in eq, and S(Z, \tau) in the call to PDESystem. Could you share the system of equations at low N, say 4? Should be output above the error.

xtalax avatar Apr 13 '22 19:04 xtalax

Here are the results at N = 4

The system of equations is:
Equation[((16.0ϕ[1](τ) + 16.0ϕ[3](τ) - 32.0ϕ[2](τ))*c[2](τ)) / ((1.0 - S(τ))^2) + IfElse.ifelse(Differential(Z)(ϕ[2](τ)) > 0, (4.0c[3](τ) - 4.0c[2](τ))*Differential(Z)(ϕ[2](τ)), (4.0c[2](τ) - 4.0c[1](τ))*Differential(Z)(ϕ[2](τ))) ~ 0.0, Differential(τ)(c[2](τ)) ~ (16.0c[1](τ) + 16.0c[3](τ) - 32.0c[2](τ)) / (0.06944444444444445((1.0 - S(τ))^2)) + (0.1790834507811324(4.0c[3](τ) - 4.0c[2](τ))*(c[2](τ)^0.5)*ϕ[2](τ)) / (0.06944444444444445 - 0.06944444444444445S(τ)), Differential(τ)(c[3](τ)) ~ (16.0c[2](τ) + 16.0c[4](τ) - 32.0c[3](τ)) / (0.06944444444444445((1.0 - S(τ))^2)) + (0.11938896718742159(4.0c[4](τ) - 4.0c[3](τ))*(c[3](τ)^0.5)*ϕ[3](τ)) / (0.06944444444444445 - 0.06944444444444445S(τ)), Differential(τ)(c[4](τ)) ~ (16.0c[3](τ) + 16.0c[5](τ) - 32.0c[4](τ)) / (0.06944444444444445((1.0 - S(τ))^2)) + (0.059694483593710795(4.0c[5](τ) - 4.0c[4](τ))*(c[4](τ)^0.5)*ϕ[4](τ)) / (0.06944444444444445 - 0.06944444444444445S(τ)), (2.0ϕ[3](τ) + 6.0ϕ[5](τ) - 8.0ϕ[4](τ))*c[5](τ) ~ 0.259102262480956S(τ) - 0.259102262480956, ϕ[5](τ) ~ 0.0, 2.0ϕ[3](τ) + 6.0ϕ[5](τ) - 8.0ϕ[4](τ) ~ 0.0, (8.0ϕ[2](τ) - 6.0ϕ[1](τ) - 2.0ϕ[3](τ)) / (1 - S(τ)) ~ -0.518204524961912(c[1](τ)^0.5) - ϕ[1](τ), 2.0c[3](τ) + 6.0c[5](τ) - 8.0c[4](τ) ~ 0.259102262480956S(τ) - 0.259102262480956, (8.0c[2](τ) - 2.0c[3](τ) - 6.0c[1](τ)) / (1 - S(τ)) ~ 0.518204524961912(c[1](τ)^0.5)*ϕ[1](τ) - 0.19398592556510205(c[1](τ)^1.5)*ϕ[1](τ)]

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Differential w.r.t. multiple variables Any[τ, Z] are not allowed.
Stacktrace:
 [1] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:171
 [2] ODESystem(deqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}}, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Any}, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, connections::Nothing, preface::Nothing, events::Vector{ModelingToolkit.SymbolicContinuousCallback}, tearing_state::Nothing, substitutions::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:120
 [3] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:174
 [4] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [6] top-level scope
   @ ~/Documents/Git/test/Scratch1.jl:142

When including Dτ(S(τ)) in eq as well as S(Z,τ) in PDESystem, an error occurs from reducing over an empty collection:

ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:301
  [2] reduce_empty(op::Function, #unused#::Type{Any})
    @ Base ./reduce.jl:311
  [3] mapreduce_empty(#unused#::typeof(identity), op::Function, T::Type)
    @ Base ./reduce.jl:345
  [4] reduce_empty(op::Base.MappingRF{typeof(identity), typeof(vcat)}, #unused#::Type{Any})
    @ Base ./reduce.jl:331
  [5] reduce_empty_iter
    @ ./reduce.jl:357 [inlined]
  [6] mapreduce_empty_iter(f::Function, op::Function, itr::Vector{Any}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:353
  [7] _mapreduce(f::typeof(identity), op::typeof(vcat), #unused#::IndexLinear, A::Vector{Any})
    @ Base ./reduce.jl:402
  [8] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::Vector{Any}, #unused#::Colon)
    @ Base ./reducedim.jl:330
  [9] #mapreduce#725
    @ ./reducedim.jl:322 [inlined]
 [10] mapreduce
    @ ./reducedim.jl:322 [inlined]
 [11] #reduce#727
    @ ./reducedim.jl:371 [inlined]
 [12] reduce(op::Function, A::Vector{Any})
    @ Base ./reducedim.jl:371
 [13] (::MethodOfLines.var"#185#189"{CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, MethodOfLines.PeriodicMap{Val{false}()}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.var"#central_ufunc#188"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}})(u::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ MethodOfLines ./none:0
 [14] iterate
    @ ./generator.jl:47 [inlined]
 [15] collect(itr::Base.Generator{Vector{Any}, MethodOfLines.var"#185#189"{CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, MethodOfLines.PeriodicMap{Val{false}()}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.var"#central_ufunc#188"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}}})
    @ Base ./array.jl:724
 [16] generate_cartesian_rules
    @ ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/generate_finite_difference_rules.jl:142 [inlined]
 [17] generate_finite_difference_rules(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, pde::Equation, derivweights::MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, pmap::MethodOfLines.PeriodicMap{Val{false}()}, indexmap::Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/generate_finite_difference_rules.jl:257
 [18] (::MethodOfLines.var"#311#313"{Equation, Term{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, Vector{Any}, MethodOfLines.PeriodicMap{Val{false}()}})(II::CartesianIndex{1})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:140
 [19] iterate
    @ ./generator.jl:47 [inlined]
 [20] _collect(c::Vector{CartesianIndex{1}}, itr::Base.Generator{Vector{CartesianIndex{1}}, MethodOfLines.var"#311#313"{Equation, Term{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, Vector{Any}, MethodOfLines.PeriodicMap{Val{false}()}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:744
 [21] collect_similar
    @ ./array.jl:653 [inlined]
 [22] map
    @ ./abstractarray.jl:2849 [inlined]
 [23] discretize_equation(pde::Equation, interior::Vector{CartesianIndex{1}}, eqvar::Term{Real, Base.ImmutableDict{DataType, Any}}, depvars::Vector{Any}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, derivweights::MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}}, indexmap::Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, boundaryvalfuncs::Vector{Any}, pmap::MethodOfLines.PeriodicMap{Val{false}()})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:138
 [24] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:98
 [25] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [26] top-level scope
    @ ~/Documents/Git/test/Scratch1.jl:143

With the updated code:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters Z τ
@variables c(..) ϕ(..) S(..)
Dz = Differential(Z)
Dzz = Differential(Z)^2
Dτ = Differential(τ)

Z_min = Z_min = τ_min = 0.0
Z_max = 1.0
τ_max = 50.0

# Parameters
iₐ = 10.
i₀ = 20.
D₁ = 1e-11
F = 96487.
R = 8.314
T = 298.15
Lₓ = 50e-6
c₀ = 1000.
Mᵥ = 1.2998e-5
Cᵣ = cᵥ = 1.

# Dimensionless Parameters
δ = iₐ*Lₓ*Cᵣ/(2*F*D₁*c₀)
tD = abs(Cᵣ)*Lₓ^2/(3600*D₁)
Da = i₀*Lₓ/(2*F*D₁*c₀)
Pe = Mᵥ*i₀*3600/(Lₓ*Cᵣ*F)
ϕ₀ = R*T/F

# PDE
eq = [
       Dz(c(Z,τ))*Dz(ϕ(Z,τ))+c(Z,τ)*(Dzz(ϕ(Z,τ)))/(1.0-S(τ))^2 ~ 0.0, 
       Dτ(c(Z,τ)) ~ (-Pe*(c(Z,τ)/c₀)^(1/2)*(-ϕ(Z,τ)/ϕ₀)*(1.0-Z)*(Dz(c(Z,τ))))/(tD*(1.0-S(τ)))+(Dzz(c(Z,τ))/(tD*(1.0-S(τ))^2)),
       Dτ(S(τ)) ~ (Pe*ϕ(Z,τ)*c(Z,τ)^(1/2))/(tD)
]
    
    
domains = [Z ∈ Interval(Z_min, Z_max),
           τ ∈ Interval(τ_min, τ_max)]

# BCs
bcs =  [
        #Initial
        c(Z,0) ~ 1.0, #Concentration @ τ=0
        S(0) ~ 0.02, #Moving Boundary @ τ=0
      
        #Boundary
        c(Z_max,τ)*(Dz(ϕ(1,τ))) ~ -δ*(1-S(τ)), #ϕ @ Z = 1
        ϕ(Z_max,τ) ~ 0.0, #ϕ @ Z = 1
        Dz(c(Z_max,τ)) ~ -δ*(1-S(τ)), #Concentration @ Z = 1
        Dz(ϕ(Z_max,τ)) ~ 0.0, #ϕ @ Z = 0
        Dz(ϕ(Z_min,τ))/(1-S(τ)) ~ -Da*(c(Z_min,τ)^(1/2))-ϕ(Z_min,τ), #ϕ BC @ Z = 0
        Dz(c(Z_min,τ))/(1-S(τ)) ~ -Da*c(Z_min,τ)^(1/2)*(-ϕ(Z_min,τ))+cᵥ*Pe*c(Z_min,τ)^(3/2)*(-ϕ(Z_min,τ)) #Concentration @ Z = 0
       ] 

@named pdesys = PDESystem(eq,bcs,domains,[Z,τ],[c(Z,τ), ϕ(Z,τ), S(Z,τ)])

N = 4

dz = (Z_max-Z_min)/N

order = 2

discretization = MOLFiniteDifference([Z => dz], τ, approx_order=order, grid_align=center_align)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

Perhaps this stems from Dτ(S(τ)) vary over τ only? The same error occurs if S(Z,τ) in the PDESystem changes to S(τ). The updated error is similar to initial error in #59, but that seemed to be solved by including the dependent variables in the PDESystem call.

BradyPlanden avatar Apr 14 '22 09:04 BradyPlanden

Ah sorry I misspoke - S(tau) is what should go in the pdesys.

I can see the problem immediately, it stems from an issue with rules matching with odd ordered derivatives, specifically the case where you have 2 odd ordered derivatives multiplying each other - defining an auxiliary variable for Dz(c(Z,τ)) with an equation dzc(Z, tau) ~ Dz(c(Z,τ)), supplying ics and bcs for this variable and replacing in the equations should be a workaround.

I will take a close look at this soon, likely will have a branch that fixes this within a week - I'll ping you when that happens. In the meantime, let me know how you get on with the workaround.

xtalax avatar Apr 14 '22 10:04 xtalax

After creating the auxiliary variable, I've encountered a similar error to #78:

ERROR: AssertionError: Variable dzc(Z, τ) does not appear in any equation, therefore cannot be solved for Stacktrace: [1] build_variable_mapping(m::Matrix{Int64}, vars::Vector{Any}, pdes::Vector{Equation}) @ MethodOfLines ~/Documents/Git/MethodOfLines.jl/src/interiormap.jl:69 [2] MethodOfLines.InteriorMap(pdes::Vector{Equation}, boundarymap::Dict{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}}}, s::MethodOfLines.DiscreteSpace{1, 4, MethodOfLines.CenterAlignedGrid}) @ MethodOfLines ~/Documents/Git/MethodOfLines.jl/src/interiormap.jl:20 [3] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}) @ MethodOfLines ~/Documents/Git/MethodOfLines.jl/src/discretization/MOL_discretization.jl:52 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}) @ MethodOfLines ~/Documents/Git/MethodOfLines.jl/src/discretization/MOL_discretization.jl:146 [5] top-level scope @ REPL[35]:2

I've also pulled the recommend branch for this issue, but ran into the same error. There appears to be nothing returned from the j = findfirst(isequal(0), cols) line in build_variable_mapping() located in interiormap.jl. This appears to stem from the creation of m from pdesys.eqs and the s = DiscreteSpace(domain, alldepvars, allindvars, discretization). Any thoughts?

Full code:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters Z τ
@variables c(..) ϕ(..) S(..) dzc(..)
Dz = Differential(Z)
Dzz = Differential(Z)^2
Dτ = Differential(τ)

Z_min = Z_min = τ_min = 0.0
Z_max = 1.0
τ_max = 50.0

# Parameters
iₐ = 10.
i₀ = 20.
D₁ = 1e-11
F = 96487.
R = 8.314
T = 298.15
Lₓ = 50e-6
c₀ = 1000.
Mᵥ = 1.2998e-5
Cᵣ = cᵥ = 1.

# Dimensionless Parameters
δ = iₐ*Lₓ*Cᵣ/(2*F*D₁*c₀)
tD = abs(Cᵣ)*Lₓ^2/(3600*D₁)
Da = i₀*Lₓ/(2*F*D₁*c₀)
Pe = Mᵥ*i₀*3600/(Lₓ*Cᵣ*F)
ϕ₀ = R*T/F

# PDE
eq = [
       dzc(Z,τ) ~ Dz(c(Z,τ)),
       dzc(Z,τ)*Dz(ϕ(Z,τ))+c(Z,τ)*(Dzz(ϕ(Z,τ)))/(1.0-S(τ))^2 ~ 0.0, 
       Dτ(c(Z,τ)) ~ (-Pe*(c(Z,τ)/c₀)^(1/2)*(-ϕ(Z,τ)/ϕ₀)*(1.0-Z)*(dzc(Z,τ)))/(tD*(1.0-S(τ)))+(Dzz(c(Z,τ))/(tD*(1.0-S(τ))^2)),
       Dτ(S(τ)) ~ (Pe*ϕ(Z,τ)*c(Z,τ)^(1/2))/(tD)
]
    
    
domains = [Z ∈ Interval(Z_min, Z_max),
           τ ∈ Interval(τ_min, τ_max)]

# BCs
bcs =  [
        #Initial
        c(Z,0) ~ 1.0, # Concentration @ τ=0
        S(0) ~ 0.02, # Moving Boundary @ τ=0
      
        #Boundary
        c(Z_max,τ)*(Dz(ϕ(1,τ))) ~ -δ*(1-S(τ)), # ϕ @ Z = 1
        ϕ(Z_max,τ) ~ 0.0, # ϕ @ Z = 1
        Dz(c(Z_max,τ)) ~ -δ*(1-S(τ)), # Concentration @ Z = 1
        Dz(ϕ(Z_max,τ)) ~ 0.0, # ϕ @ Z = 0
        Dz(ϕ(Z_min,τ))/(1-S(τ)) ~ -Da*(c(Z_min,τ)^(1/2))-ϕ(Z_min,τ), # ϕ BC @ Z = 0
        Dz(c(Z_min,τ))/(1-S(τ)) ~ -Da*c(Z_min,τ)^(1/2)*(-ϕ(Z_min,τ))+cᵥ*Pe*c(Z_min,τ)^(3/2)*(-ϕ(Z_min,τ)), # Concentration @ Z = 0
        dzc(Z_max,τ) ~ -δ*(1-S(τ)),
        dzc(Z_min,τ)/(1-S(τ)) ~ -Da*c(Z_min,τ)^(1/2)*(-ϕ(Z_min,τ))+cᵥ*Pe*c(Z_min,τ)^(3/2)*(-ϕ(Z_min,τ)) # Concentration @ Z = 0
       ] 

@named pdesys = PDESystem(eq,bcs,domains,[Z,τ],[c(Z,τ), ϕ(Z,τ), S(τ), dzc(Z,τ)])

N = 4

dz = (Z_max-Z_min)/N

order = 2

discretization = MOLFiniteDifference([Z => dz], τ, approx_order=order, grid_align=center_align)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

BradyPlanden avatar Apr 16 '22 12:04 BradyPlanden

I get a different error on that branch - I'll take a closer look anyway

xtalax avatar Apr 25 '22 15:04 xtalax

I've tracked the issue to this equation Dτ(S(τ)) ~ (Pe*ϕ(Z,τ)*c(Z,τ)^(1/2))/(tD). I don't know what to do about this since the lhs has no Z dependence, but the rhs does, can you advise on how this would be discretized?

xtalax avatar Apr 25 '22 17:04 xtalax