From f4139aabcd1e8ba2c3a047906f2700bad41f0f1b Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Wed, 18 Oct 2023 15:40:01 +0100 Subject: [PATCH 1/6] a --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c4b4b9719..2db1a2f8d 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,4 @@ dev/* # Build artifacts for creating documentation generated by the Documenter package docs/build/ docs/site/ +.DS_Store \ No newline at end of file From 422e0a0daf3489cc092946e8c539b0d8b6df9a68 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 19 Oct 2023 17:21:57 +0100 Subject: [PATCH 2/6] add callbacks --- src/MethodOfLines.jl | 4 ++- .../differential_discretizer.jl | 4 ++- .../generate_finite_difference_rules.jl | 4 ++- .../schemes/callbacks/callback_rules.jl | 5 +++ src/interface/MOLFiniteDifference.jl | 5 +-- src/interface/callbacks.jl | 13 +++++++ test/components/callbacks.jl | 35 +++++++++++++++++++ 7 files changed, 65 insertions(+), 5 deletions(-) create mode 100644 src/discretization/schemes/callbacks/callback_rules.jl create mode 100644 src/interface/callbacks.jl create mode 100644 test/components/callbacks.jl diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 7a818230e..5bd8e6107 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -50,6 +50,7 @@ import Base.ndims # Interface include("interface/grid_types.jl") include("interface/scheme_types.jl") +include("interface/callbacks.jl") include("interface/disc_strategy_types.jl") include("interface/MOLFiniteDifference.jl") @@ -72,6 +73,7 @@ include("discretization/schemes/upwind_difference/upwind_diff_weights.jl") include("discretization/schemes/half_offset_weights.jl") include("discretization/schemes/extrapolation_weights.jl") include("discretization/differential_discretizer.jl") +include("discretization/schemes/callbacks/callback_rules.jl") # System Parsing include("system_parsing/pde_system_transformation.jl") @@ -104,6 +106,6 @@ include("precompile.jl") # Export export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete, chebyspace -export UpwindScheme, WENOScheme, FunctionalScheme +export UpwindScheme, WENOScheme, FunctionalScheme, MOLDiscCallback end diff --git a/src/discretization/differential_discretizer.jl b/src/discretization/differential_discretizer.jl index 5a6250c76..d23615d0e 100644 --- a/src/discretization/differential_discretizer.jl +++ b/src/discretization/differential_discretizer.jl @@ -8,6 +8,7 @@ struct DifferentialDiscretizer{T, D1, S} <: AbstractDifferentialDiscretizer interpmap::Dict{Num,DerivativeOperator} orders::Dict{Num,Vector{Int}} boundary::Dict{Num,DerivativeOperator} + callbacks end function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, discretization::MOLFiniteDifference, orders) @@ -15,6 +16,7 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di bcs = pdesys.bcs isa Vector ? pdesys.bcs : [pdesys.bcs] approx_order = discretization.approx_order advection_scheme = discretization.advection_scheme + callbacks = discretization.callbacks upwind_order = advection_scheme isa UpwindScheme ? advection_scheme.order : 1 differentialmap = Array{Pair{Num,DerivativeOperator},1}() @@ -63,5 +65,5 @@ function PDEBase.construct_differential_discretizer(pdesys, s::DiscreteSpace, di push!(boundary, x => BoundaryInterpolatorExtrapolator(max(6, approx_order), dx)) end - return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary)) + return DifferentialDiscretizer{eltype(orders),typeof(Dict(differentialmap)),typeof(advection_scheme)}(approx_order, advection_scheme, Dict(differentialmap), (Dict(nonlinlap_inner), Dict(nonlinlap_outer)), (Dict(windpos), Dict(windneg)), Dict(interp), Dict(orders), Dict(boundary), callbacks) end diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 8d2ad1b5d..e5eacab33 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -80,6 +80,8 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, integration_rules = [] end + cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) - return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) + return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules, cb_rules) end diff --git a/src/discretization/schemes/callbacks/callback_rules.jl b/src/discretization/schemes/callbacks/callback_rules.jl new file mode 100644 index 000000000..cc88217e5 --- /dev/null +++ b/src/discretization/schemes/callbacks/callback_rules.jl @@ -0,0 +1,5 @@ +function generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + cbs = derivweights.callbacks + cb_rules = [cb.sym => cb(s) for cb in cbs] + return cb_rules +end \ No newline at end of file diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl index a4a7d7173..968e105e2 100644 --- a/src/interface/MOLFiniteDifference.jl +++ b/src/interface/MOLFiniteDifference.jl @@ -33,11 +33,12 @@ struct MOLFiniteDifference{G,D} <: AbstractEquationSystemDiscretization use_ODAE::Bool disc_strategy::D useIR::Bool + callbacks kwargs end # Constructors. If no order is specified, both upwind and centered differences will be 2nd order -function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, kwargs...) +function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_scheme = UpwindScheme(), grid_align=CenterAlignedGrid(), discretization_strategy = ScalarizedDiscretization(), upwind_order = nothing, should_transform = true, use_ODAE = false, useIR = true, callbacks = [], kwargs...) if upwind_order !== nothing @warn "`upwind_order` no longer does anything, and will be removed in a future release. See the docs for the current interface." end @@ -50,7 +51,7 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche dxs = dxs isa Dict ? dxs : Dict(dxs) - return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, kwargs) + return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, useIR, callbacks, kwargs) end PDEBase.get_time(disc::MOLFiniteDifference) = disc.time diff --git a/src/interface/callbacks.jl b/src/interface/callbacks.jl new file mode 100644 index 000000000..2683e1a4f --- /dev/null +++ b/src/interface/callbacks.jl @@ -0,0 +1,13 @@ +struct MOLDiscCallback + f::Function + p + sym +end + +function MOLDiscCallback(s, disc_ps) + symh = hash(f) + sym = unwrap(first(@parameters(Symbol("MOLDiscCallback_$symh")))) + return MOLDiscCallback(f, disc_ps, sym) +end + +(M::MOLDiscCallback)(s) = M.f(s, M.p) \ No newline at end of file diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl new file mode 100644 index 000000000..1a8721217 --- /dev/null +++ b/test/components/callbacks.jl @@ -0,0 +1,35 @@ +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils + +@testset "Discrete callback" begin + @parameters x, t + @variables u(..) + Dt = Differential(t) + Dxx = Differential(x)^2 + t_min = 0.0 + t_max = 2.0 + x_min = 0.0 + x_max = 2.0 + + a = 0.1*pi + + eq = Dt(u(t, x)) ~ a*Dxx(u(t, x)) + + bcs = [u(t_min, x) ~ sinpi(x), u(t, x_min) ~ sinpi(t+x), u(t, x_max) ~ sinpi(t+x)] + + domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] + + pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + + cb = MOLDiscCallback((s, p) -> prod(p), [0.1, pi]) + a_cb = cb.sym + cbeq = Dt(u(t, x)) ~ a_cb*Dxx(u(t, x)) + + cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)], callbacks = [cb]) + + disc = MOLFiniteDifference([x => 0.1], t; approx_order=2) + + prob = discretize(pdesys, disc) + cbprob = discretize(cbpdesys, disc) + + @test solve(prob, Tsit5(), saveat=0.1) ≈ solve(cbprob, Tsit5(), saveat=0.1) +end \ No newline at end of file From ac10b6af34e78726fc9577fd0d8f4bc77ff19c20 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Thu, 19 Oct 2023 17:38:43 +0100 Subject: [PATCH 3/6] runtests --- src/discretization/generate_finite_difference_rules.jl | 2 +- test/runtests.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index e5eacab33..ac11911ec 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -83,5 +83,5 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, cb_rules = generate_cb_rules(II, s, depvars, derivweights, bmap, indexmap, terms) integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) - return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules, cb_rules) + return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) end diff --git a/test/runtests.jl b/test/runtests.jl index 257e158fd..263692a8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -71,6 +71,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") @time @safetestset "ODEFunction" begin include("components/ODEFunction_test.jl") end + @time @safetestset "Discrete Callbacks" begin + include("components/callbacks.jl") + end #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end end From ccc47f9fd53aec7c3b02b10ccd38818da239f6ba Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 18:47:46 +0000 Subject: [PATCH 4/6] callbacks passing --- src/interface/callbacks.jl | 2 +- src/interface/solution/timeindep.jl | 2 +- src/system_parsing/pde_system_transformation.jl | 3 +++ test/components/callbacks.jl | 13 ++++++++----- 4 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/interface/callbacks.jl b/src/interface/callbacks.jl index 2683e1a4f..5c5fc427e 100644 --- a/src/interface/callbacks.jl +++ b/src/interface/callbacks.jl @@ -4,7 +4,7 @@ struct MOLDiscCallback sym end -function MOLDiscCallback(s, disc_ps) +function MOLDiscCallback(f, disc_ps) symh = hash(f) sym = unwrap(first(@parameters(Symbol("MOLDiscCallback_$symh")))) return MOLDiscCallback(f, disc_ps, sym) diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl index 8e3374952..b1239c845 100644 --- a/src/interface/solution/timeindep.jl +++ b/src/interface/solution/timeindep.jl @@ -11,7 +11,7 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada umap = mapreduce(vcat, dvs) do u let discu = discretespace.discvars[u] solu = map(CartesianIndices(discu)) do I - i = sym_to_index(discu[I], odesys.states) + i = sym_to_index(discu[I], get_states(odesys)) # Handle Observed if i !== nothing sol.u[i] diff --git a/src/system_parsing/pde_system_transformation.jl b/src/system_parsing/pde_system_transformation.jl index 681b00866..ca722b151 100644 --- a/src/system_parsing/pde_system_transformation.jl +++ b/src/system_parsing/pde_system_transformation.jl @@ -35,6 +35,9 @@ function PDEBase.transform_pde_system!(v::PDEBase.VariableMap, boundarymap, sys: end function PDEBase.should_transform(pdesys::PDESystem, disc::MOLFiniteDifference, boundarymap) + if !disc.should_transform + return false + end if has_interfaces(boundarymap) @warn "The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature." return false diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl index 1a8721217..a1d1236f9 100644 --- a/test/components/callbacks.jl +++ b/test/components/callbacks.jl @@ -1,4 +1,4 @@ -using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, OrdinaryDiffEq @testset "Discrete callback" begin @parameters x, t @@ -18,18 +18,21 @@ using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)] - pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) + @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)]) cb = MOLDiscCallback((s, p) -> prod(p), [0.1, pi]) a_cb = cb.sym cbeq = Dt(u(t, x)) ~ a_cb*Dxx(u(t, x)) - cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)], callbacks = [cb]) + @named cbpdesys = PDESystem(cbeq, bcs, domains, [t, x], [u(t, x)]) - disc = MOLFiniteDifference([x => 0.1], t; approx_order=2) + disc = MOLFiniteDifference([x => 0.1], t; approx_order=2, callbacks = [cb]) prob = discretize(pdesys, disc) cbprob = discretize(cbpdesys, disc) - @test solve(prob, Tsit5(), saveat=0.1) ≈ solve(cbprob, Tsit5(), saveat=0.1) + sol1 = solve(prob, Tsit5(), saveat=0.1) + sol2 = solve(cbprob, Tsit5(), saveat=0.1) + + @test sol1[u(t, x)] ≈ sol2[u(t, x)] end \ No newline at end of file From 135c86a51f6b2a75aa01f9129d81c8bfa94c0c83 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 19:55:15 +0000 Subject: [PATCH 5/6] test --- test/components/callbacks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/components/callbacks.jl b/test/components/callbacks.jl index a1d1236f9..e5a6b89ea 100644 --- a/test/components/callbacks.jl +++ b/test/components/callbacks.jl @@ -5,7 +5,7 @@ using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils @variables u(..) Dt = Differential(t) Dxx = Differential(x)^2 - t_min = 0.0 + t_min = 0.0 t_max = 2.0 x_min = 0.0 x_max = 2.0 From 7966aa2e15b33730135e6ed2a32c26da64ff852e Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 14 Nov 2023 19:55:59 +0000 Subject: [PATCH 6/6] update compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e897ccd0f..c1b21a2a8 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ Interpolations = "0.14" Latexify = "0.15, 0.16" ModelingToolkit = "8" OrdinaryDiffEq = "6" -PDEBase = "0.1.4" +PDEBase = "0.1.6" PrecompileTools = "1" RuntimeGeneratedFunctions = "0.5.9" SafeTestsets = "0.0.1"