I would expect discretization to work correctly when units are involved.
using ModelingToolkit
using ModelingToolkit: t, D
using MethodOfLines, DomainSets, DynamicQuantities, OrdinaryDiffEqDefault
@parameters x [unit = u"m"]
@parameters D_diff [description = "Diffusion coefficient", unit = u"m^2/s"]
@variables u(..) [unit = u"kg/m^3"]
@constants begin
u_ref = 1.0, [description = "Reference concentration", unit = u"kg/m^3"]
x_ref = 1.0, [description = "Reference length", unit = u"m"]
end
Dx = Differential(x)
Dxx = Dx ∘ Dx # Second derivative operator
eq = [D(u(t, x)) ~ D_diff * Dxx(u(t, x))] # Continuous PDE form
# Define spatial domain and boundary/initial conditions
domains = [t ∈ Interval(0.0, 1.0),
x ∈ Interval(0.0, 1.0)]
bcs = [u(0, x) ~ u_ref * sin(π * x / x_ref), # Initial condition
u(t, 0) ~ 0 * u_ref, # Left boundary
u(t, 1) ~ 0 * u_ref] # Right boundary
@named sys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [D_diff => 0.1])
dx = 0.05
discretization = MOLFiniteDifference([x => dx], t)
prob = MethodOfLines.discretize(sys, discretization)
sol = solve(prob)
┌ Warning: in eq. #1left, in sum Differential(t)((u(t))[2]) - D_diff*(399.99999999999994(u(t))[1] - 799.9999999999999(u(t))[2] + 399.99999999999994(u(t))[3]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #2left, in sum Differential(t)((u(t))[3]) - D_diff*(399.99999999999994(u(t))[2] - 799.9999999999999(u(t))[3] + 399.99999999999994(u(t))[4]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #3left, in sum Differential(t)((u(t))[4]) - D_diff*(399.99999999999994(u(t))[3] - 799.9999999999999(u(t))[4] + 399.99999999999994(u(t))[5]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #4left, in sum Differential(t)((u(t))[5]) - D_diff*(399.99999999999994(u(t))[4] - 799.9999999999999(u(t))[5] + 399.99999999999994(u(t))[6]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #5left, in sum Differential(t)((u(t))[6]) - D_diff*(399.99999999999994(u(t))[5] - 799.9999999999999(u(t))[6] + 399.99999999999994(u(t))[7]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #6left, in sum Differential(t)((u(t))[7]) - D_diff*(399.99999999999994(u(t))[6] - 799.9999999999999(u(t))[7] + 399.99999999999994(u(t))[8]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #7left, in sum Differential(t)((u(t))[8]) - D_diff*(399.99999999999994(u(t))[7] - 799.9999999999999(u(t))[8] + 399.99999999999994(u(t))[9]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #8left, in sum Differential(t)((u(t))[9]) - D_diff*(399.99999999999994(u(t))[10] + 399.99999999999994(u(t))[8] - 799.9999999999999(u(t))[9]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #9left, in sum Differential(t)((u(t))[10]) - D_diff*(-799.9999999999999(u(t))[10] + 399.99999999999994(u(t))[11] + 399.99999999999994(u(t))[9]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #10left, in sum Differential(t)((u(t))[11]) - D_diff*(399.99999999999994(u(t))[10] - 799.9999999999999(u(t))[11] + 399.99999999999994(u(t))[12]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #11left, in sum Differential(t)((u(t))[12]) - D_diff*(399.99999999999994(u(t))[11] - 799.9999999999999(u(t))[12] + 399.99999999999994(u(t))[13]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #12left, in sum Differential(t)((u(t))[13]) - D_diff*(399.99999999999994(u(t))[12] - 799.9999999999999(u(t))[13] + 399.99999999999994(u(t))[14]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #13left, in sum Differential(t)((u(t))[14]) - D_diff*(399.99999999999994(u(t))[13] - 799.9999999999999(u(t))[14] + 399.99999999999994(u(t))[15]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #14left, in sum Differential(t)((u(t))[15]) - D_diff*(399.99999999999994(u(t))[14] - 799.9999999999999(u(t))[15] + 399.99999999999994(u(t))[16]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #15left, in sum Differential(t)((u(t))[16]) - D_diff*(399.99999999999994(u(t))[15] - 799.9999999999999(u(t))[16] + 399.99999999999994(u(t))[17]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #16left, in sum Differential(t)((u(t))[17]) - D_diff*(399.99999999999994(u(t))[16] - 799.9999999999999(u(t))[17] + 399.99999999999994(u(t))[18]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #17left, in sum Differential(t)((u(t))[18]) - D_diff*(399.99999999999994(u(t))[17] - 799.9999999999999(u(t))[18] + 399.99999999999994(u(t))[19]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #18left, in sum Differential(t)((u(t))[19]) - D_diff*(399.99999999999994(u(t))[18] - 799.9999999999999(u(t))[19] + 399.99999999999994(u(t))[20]), units [1.0 m² s⁻¹, 1.0 s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
┌ Warning: in eq. #19left, in sum Differential(t)((u(t))[20]) - D_diff*(399.99999999999994(u(t))[19] - 799.9999999999999(u(t))[20] + 399.99999999999994(u(t))[21]), units [1.0 s⁻¹, 1.0 m² s⁻¹] do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:181
The system of equations is:
Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.
ERROR: ModelingToolkit.ValidationError("Some equations had invalid units. See warnings for details.")
Stacktrace:
[1] check_units
@ ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/unit_check.jl:329 [inlined]
[2] System(tag::UInt64, eqs::Vector{…}, noise_eqs::Nothing, jumps::Vector{…}, constraints::Vector{…}, costs::Vector{…}, consolidate::typeof(ModelingToolkit.default_consolidate), unknowns::Vector{…}, ps::Vector{…}, brownians::Vector{…}, iv::SymbolicUtils.BasicSymbolic{…}, observed::Vector{…}, parameter_dependencies::Vector{…}, var_to_name::Dict{…}, name::Symbol, description::String, defaults::Dict{…}, guesses::Dict{…}, systems::Vector{…}, initialization_eqs::Vector{…}, continuous_events::Vector{…}, discrete_events::Vector{…}, connector_type::Nothing, assertions::Dict{…}, metadata::Base.ImmutableDict{…}, gui_metadata::Nothing, is_dde::Bool, tstops::Vector{…}, inputs::OrderedCollections.OrderedSet{…}, outputs::OrderedCollections.OrderedSet{…}, tearing_state::Nothing, namespacing::Bool, complete::Bool, index_cache::Nothing, ignored_connections::Nothing, preface::Vector{…}, parent::Nothing, initializesystem::Nothing, is_initializesystem::Bool, is_discrete::Bool, isscheduled::Bool, schedule::Nothing; checks::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/system.jl:296
[3] System
@ ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/system.jl:263 [inlined]
[4] System(eqs::Vector{…}, iv::Num, dvs::Vector{…}, ps::Vector{…}, brownians::Vector{…}; constraints::Vector{…}, noise_eqs::Nothing, jumps::Vector{…}, costs::Vector{…}, consolidate::Function, observed::Vector{…}, parameter_dependencies::Vector{…}, defaults::Dict{…}, guesses::Dict{…}, systems::Vector{…}, initialization_eqs::Vector{…}, continuous_events::Vector{…}, discrete_events::Vector{…}, connector_type::Nothing, assertions::Dict{…}, metadata::Vector{…}, gui_metadata::Nothing, is_dde::Nothing, tstops::Vector{…}, inputs::OrderedCollections.OrderedSet{…}, outputs::OrderedCollections.OrderedSet{…}, tearing_state::Nothing, ignored_connections::Nothing, parent::Nothing, description::String, name::Symbol, discover_from_metadata::Bool, initializesystem::Nothing, is_initializesystem::Bool, is_discrete::Bool, preface::Vector{…}, checks::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/JvjlW/src/systems/system.jl:452
[5] generate_system(disc_state::PDEBase.EquationState, s::MethodOfLines.DiscreteSpace{…}, u0::Vector{…}, tspan::Tuple{…}, metadata::MethodOfLines.MOLMetadata{…}, disc::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/CpWUp/src/discretization_state.jl:40
[6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/CpWUp/src/symbolic_discretize.jl:90
[7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase ~/.julia/packages/PDEBase/CpWUp/src/discretization_state.jl:57
[8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/CpWUp/src/discretization_state.jl:54
[9] top-level scope
@...
Some type information was truncated. Use `show(err)` to see complete types.
Describe the bug 🐞
A duplicate of #287 , but with an updated example
Expected behavior
I would expect discretization to work correctly when units are involved.
Minimal Reproducible Example 👇
Error & Stacktrace⚠️
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()