Comments (8)
ChrisRackauckas
commented on June 6, 2024
Care you share a runnable example?
from methodoflines.jl.
SimonEnsemble
commented on June 6, 2024
sorry @ChrisRackauckas, will do for now on.
here is a minimal example. see the four different WAYs. all are the same mathematically, but two yield the error above.
v4.25.0 in DiffEqOperators.
using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, DataFrames, CSV, Statistics, PlutoUI, DifferentialEquations
function toy_pde()
@parameters x t
@variables c(..)
∂t = Differential(t)
∂x = Differential(x)
∂²x = Differential(x) ^ 2
D₀ = 1.5
α = 0.15
χ = 1.2
R = 0.1
cₑ = 2.0
ℓ = 1.0
Δx = 0.01
# WAY 1: MethodError
# D = D₀ / (1.0 + exp(α * (c(x, t) - χ)))
# diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
# WAY 2: no error. it's the D₀ then?
# D = 1.0 / (1.0 + exp(α * (c(x, t) - χ)))
# diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
# WAY 3: no error.
# diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))
# WAY 4: error.
diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))
bcs = [
# initial condition
c(x, 0) ~ 0.0,
# Robin BC
∂x(c(0.0, t)) / (1 + exp(α * (c(0.0, t) - χ))) * R * D₀ + cₑ - c(0.0, t) ~ 0.0,
# no flux BC
∂x(c(ℓ, t)) ~ 0.0]
# define space-time plane
domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]
# put it all together into a PDE system
@named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]);
# discretize space
discretization = MOLFiniteDifference([x=>Δx], t)
# convert the PDE into a system of ODEs via method of lines
prob = discretize(pdesys, discretization)
# compute number of spatial discretization points, boundaries included.
Nₓ = length(prob.u0) + 1 # change to +2 for Dirichlet
@assert (length(0:Δx:ℓ) == Nₓ)
# solve system of ODEs in time.
sol = solve(prob, saveat=0.01)
# convert solution to data frame. solution at boundaries must be inferred.
return DataFrame(sol)
end
toy_pde()
from methodoflines.jl.
ChrisRackauckas
commented on June 6, 2024
Thanks a ton. That'll be good debugging fodder.
from methodoflines.jl.
SimonEnsemble
commented on June 6, 2024
@ChrisRackauckas here is another example with a concentration-dependent diffusion coefficient that fails with an error. it works until discretize is called.
using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, CairoMakie,ColorSchemes, DataFrames, CSV, Statistics, Optim, Polynomials, PlutoUI, DifferentialEquations, SpecialFunctions
function toy_pde2()
@parameters x t
@variables c(..)
∂t = Differential(t)
∂x = Differential(x)
∂²x = Differential(x) ^ 2
D₀ = 1.2
c_max = 0.15
c_e = 0.1
η = 0.1
cₑ = 2.0
ℓ = 1.0
Δx = 0.01
# D = D(c)
θ = c(x, t) / c_max
β = sqrt(1 - 4 * θ * (1 - θ) * (1 - η))
ϵ = (2 * θ + β - 1) / (2 * η * (1 - θ))
D = D₀ * (1 + ϵ) ^ 3 / (1 - θ) / (1 + η * ϵ) ^ 4 * (1 + 2 * (1 - β) / β)
diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
bcs = [
# initial condition
c(x, 0) ~ 0.0,
# Dirichelt BC
c(0.0, t) ~ c_e,
# no flux BC
∂x(c(ℓ, t)) ~ 0.0]
# define space-time plane
domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]
# put it all together into a PDE system
@named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]);
# discretize space
discretization = MOLFiniteDifference([x=>Δx], t)
# convert the PDE into a system of ODEs via method of lines
prob = discretize(pdesys, discretization)
end
toy_pde2()the error is:
MethodError: no method matching collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Div{Real, SymbolicUtils.Mul{Real, Float64, Dict{Any, Number}, Nothing}, SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, Nothing})
Closest candidates are:
collect_ivs_from_nested_differential!(::Any, !Matched::SymbolicUtils.Term) at ~/.julia/packages/ModelingToolkit/r5ZaU/src/utils.jl:160
collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Term{Real, Nothing})@utils.jl:164
collect_differentials(::Vector{Symbolics.Equation})@utils.jl:138
check_equations(::Vector{Symbolics.Equation}, ::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}})@utils.jl:151
(::ModelingToolkit.var"#ODESystem#92#93")(::Bool, ::Type{ModelingToolkit.ODESystem}, ::Vector{Symbolics.Equation}, ::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, ::Vector{SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}}}, ::Vector{Any}, ::Dict{Any, Any}, ::Vector{Any}, ::Vector{Symbolics.Equation}, ::Base.RefValue{Vector{Symbolics.Num}}, ::Base.RefValue{Any}, ::Base.RefValue{Any}, ::Base.RefValue{Matrix{Symbolics.Num}}, ::Base.RefValue{Matrix{Symbolics.Num}}, ::Symbol, ::Vector{ModelingToolkit.ODESystem}, ::Dict{Any, Any}, ::Nothing, ::Nothing, ::Nothing, ::Vector{ModelingToolkit.SymbolicContinuousCallback})@odesystem.jl:100
var"#ODESystem#94"(::Vector{Symbolics.Num}, ::Vector{Symbolics.Equation}, ::Vector{ModelingToolkit.ODESystem}, ::Symbol, ::Dict{Any, Any}, ::Dict{Any, Any}, ::Dict{Symbolics.Num, Float64}, ::Nothing, ::Nothing, ::Nothing, ::Bool, ::Type{ModelingToolkit.ODESystem}, ::Vector{Symbolics.Equation}, ::Symbolics.Num, ::Vector{Symbolics.Num}, ::Vector{Symbolics.Num})@odesystem.jl:151
symbolic_discretize(::ModelingToolkit.PDESystem, ::DiffEqOperators.MOLFiniteDifference{Vector{Pair{Symbolics.Num, Float64}}, Symbolics.Num})@MOL_discretization.jl:400
discretize(::ModelingToolkit.PDESystem, ::DiffEqOperators.MOLFiniteDifference{Vector{Pair{Symbolics.Num, Float64}}, Symbolics.Num})@MOL_discretization.jl:406
toy_pde2()@Other: 45
top-level scope@Local: 1[inlined]
from methodoflines.jl.
SimonEnsemble
commented on June 6, 2024
we've isolated the problem as due to division:
ϵ = (2 * θ + β - 1) # / (2 * η * (1 - θ))
D = D₀ * (1 + ϵ) ^ 3# / (1 - θ) / (1 + η * ϵ) ^ 4 * (1 + 2 * (1 - β) / β)prevents the error.
from methodoflines.jl.
SimonEnsemble
commented on June 6, 2024
the key was to add expand_derivatives().
diff_eq = expand_derivatives(∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))))makes the error disappear. found this in:
SciML/ModelingToolkit.jl#1197
from methodoflines.jl.
xtalax
commented on June 6, 2024
This is due to limitations in @rule, in particular the lack of equivalance between *(a^-1, u(t, x), b and /(*(b, u(t, x), a) in the matching. a future PR will fix this, but this gets complicated very quickly as you have to match each possible set of nested functions individually
A workaround would be to expand / in to a flat equivalent * form before matching, but this doesn't avoid the issue of when some u is nested inside a function like exp
from methodoflines.jl.
xtalax
commented on June 6, 2024
@shashi What should I do about this?
from methodoflines.jl.
Related Issues (20)
- Suppress output in discretization HOT 1
- Unstable steady-state solutions for smaller grid sizes in heat equation example HOT 1
- Fail at discretization HOT 1
- Discretization failure: MOLFiniteDifference-created variables are not included in variable map when creating ODEProblem HOT 2
- Error with Solving Differential Equations HOT 2
- Parallel computing HOT 2
- Spherical laplacian example seems to give moderate error between exact and numerical solution
- Issue with forwarddiff
- FAQ
- @parameters work unexpectedly HOT 3
- Units don't seem to work correctly HOT 4
- Auxiliary variable case fails
- PIDE Tutorial HOT 1
- 'update_varmap!' not being used HOT 3
- ERROR: type Array has no field lhs for Schrodinger equation HOT 2
- Can a parameter be used as a domain boundary? HOT 4
- Ifelse condition on a BC HOT 1
- Discretization failed;ArgumentError: Differential w.r.t. multiple variables are not allowed.
- Jet Bundle Power Series Discretization HOT 1
- Misleading warning in presence of auxiliary variables
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from methodoflines.jl.