diff --git a/src/interface/solution/common.jl b/src/interface/solution/common.jl index 75762e921..96e16abd2 100644 --- a/src/interface/solution/common.jl +++ b/src/interface/solution/common.jl @@ -13,7 +13,7 @@ function (sol::SciMLBase.PDESolution{T,N,S,D})(args...; if dv === nothing @assert length(args) == length(sol.ivs) "Not enough arguments for the number of independent variables including time where appropriate, got $(length(args)) expected $(length(sol.ivs))." return map(sol.dvs) do dv - arg_ivs = arguments(dv.val) + arg_ivs = arguments(safe_unwrap(dv)) is = map(arg_ivs) do arg_iv i = findfirst(isequal(arg_iv), sol.ivs) @assert i !== nothing "Independent variable $(arg_iv) in dependent variable $(dv) not found in the solution." @@ -37,6 +37,8 @@ Base.@propagate_inbounds function Base.getindex(A::SciMLBase.PDESolution{T,N,S,D idv = sym_to_index(sym, A.dvs) if idv !== nothing dv = A.dvs[idv] + elseif any(isequal(safe_unwrap(sym)), safe_unwrap.(collect(values(A.disc_data.discretespace.vars.replaced_vars)))) + dv = sym end if SciMLBase.issymbollike(sym) && iv !== nothing && isequal(sym, iv) A.ivdomain[iiv] diff --git a/src/interface/solution/solution_utils.jl b/src/interface/solution/solution_utils.jl index 99d877be0..c1a641606 100644 --- a/src/interface/solution/solution_utils.jl +++ b/src/interface/solution/solution_utils.jl @@ -1,7 +1,7 @@ # TODO: Check if the grid is uniform and use the supported higher order interpolations -function build_interpolation(umap, ivs, ivgrid, sol, pdesys) - return Dict(map(collect(keys(umap))) do k +function build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys) + return Dict(map(Num.(dvs)) do k args = arguments(k.val) nodes = (map(args) do arg i = findfirst(isequal(arg), ivs) diff --git a/src/interface/solution/timedep.jl b/src/interface/solution/timedep.jl index f0c731ed6..f52812f32 100644 --- a/src/interface/solution/timedep.jl +++ b/src/interface/solution/timedep.jl @@ -14,8 +14,9 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, else states(odesys) end + dvs = discretespace.ū # Reshape the solution to flat arrays, faster to do this eagerly. - umap = Dict(map(discretespace.ū) do u + umap = mapreduce(vcat, dvs) do u let discu = discretespace.discvars[u] solu = map(CartesianIndices(discu)) do I i = sym_to_index(discu[I], solved_states) @@ -38,14 +39,20 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, out[I, :] .= solu[I] end else - @assert false "The time variable must be the first or last argument of the dependent variable $u." + error("The time variable must be the first or last argument of the dependent variable $u.") end - Num(u) => out + # Deal with any replaced variables + ureplaced = get(discretespace.vars.replaced_vars, u, nothing) + if isnothing(ureplaced) + [Num(u) => out] + else + [Num(u) => out, ureplaced => out] + end end - end) + end |> Dict # Build Interpolations - interp = build_interpolation(umap, ivs, ivgrid, sol, pdesys) + interp = build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys) return SciMLBase.PDETimeSeriesSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata), typeof(sol),typeof(sol.errors),typeof(sol.t),typeof(ivgrid), diff --git a/src/interface/solution/timeindep.jl b/src/interface/solution/timeindep.jl index 3fce94195..691c8b873 100644 --- a/src/interface/solution/timeindep.jl +++ b/src/interface/solution/timeindep.jl @@ -6,8 +6,9 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada # Extract axies ivs = [discretespace.x̄...] ivgrid = ((discretespace.grid[x] for x in ivs)...,) + dvs = discretespace.ū # Reshape the solution to flat arrays - umap = Dict(map(discretespace.ū) do u + 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) @@ -22,11 +23,17 @@ function SciMLBase.PDENoTimeSolution(sol::SciMLBase.NonlinearSolution{T}, metada for I in CartesianIndices(discu) out[I] = solu[I] end - Num(u) => out + # Deal with any replaced variables + ureplaced = get(discretespace.vars.replaced_vars, u, nothing) + if isnothing(ureplaced) + [Num(u) => out] + else + [Num(u) => out, ureplaced => out] + end end - end) + end |> Dict # Build Interpolations - interp = build_interpolation(umap, ivs, ivgrid, sol, pdesys) + interp = build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys) return SciMLBase.PDENoTimeSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata), typeof(sol),typeof(ivgrid),typeof(ivs),typeof(pdesys.dvs),typeof(sol.prob),typeof(sol.alg), diff --git a/test/pde_systems/MOLtest2.jl b/test/pde_systems/MOLtest2.jl index cef11ec27..0dc1892df 100644 --- a/test/pde_systems/MOLtest2.jl +++ b/test/pde_systems/MOLtest2.jl @@ -38,11 +38,12 @@ using Test # define space-time plane domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)] - @testset "Test 01: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin + @test_broken begin #@testset "Test 01: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin D = D₀ / (1.0 + exp(α * (c(x, t) - χ))) diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end @testset "Test 02: ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))" begin @@ -50,60 +51,70 @@ using Test diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end @testset "Test 03: ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0 / D₀ + exp(α * (c(x, t) - χ)) / D₀) * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end - @testset "Test 04: ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))" begin + @test_broken begin #@testset "Test 04: ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end @testset "Test 05: ∂t(c(x, t)) ~ ∂x(1/x * ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ ∂x(1 / x * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end - @testset "Test 06: ∂t(c(x, t)) ~ ∂x(x*∂x(c(x, t)))/c(x,t)" begin + @test_broken begin #@testset "Test 06: ∂t(c(x, t)) ~ ∂x(x*∂x(c(x, t)))/c(x,t)" begin diff_eq = ∂t(c(x, t)) ~ ∂x(x * ∂x(c(x, t))) / c(x, t) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end @testset "Test 07: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)) ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)) * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end - @testset "Test 08: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))" begin + @test_broken begin #@testset "Test 08: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end - @testset "Test 09: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin - diff_eq = c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) / (1 + c(x, t)) + @test_broken begin #@testset "Test 09: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin + diff_eq = c(x, t) * ∂x(c(x, t) * ∂x(c(x, t))) / (1 + c(x, t)) ~ 0 @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end - @testset "Test 10: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin - diff_eq = c(x, t) * ∂x(c(x, t)^(-1) * ∂x(c(x, t))) + @test_broken begin #@testset "Test 10: ∂t(c(x, t)) ~ c(x, t) * ∂x(c(x,t) * ∂x(c(x, t)))/(1+c(x,t))" begin + diff_eq = c(x, t) * ∂x(c(x, t)^(-1) * ∂x(c(x, t))) ~ 0 @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end @testset "Test 11: ∂t(c(x, t)) ~ ∂x(1/(1+c(x,t)^2) ∂x(c(x, t)))" begin diff_eq = ∂t(c(x, t)) ~ ∂x(1 / (1 + c(x, t)^2) * ∂x(c(x, t))) @named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]) discretization = MOLFiniteDifference([x => Δx], t) + prob = discretize(pdesys, discretization) end end @@ -361,3 +372,53 @@ end solu = sol[u(t, x)] solv = sol[v(t)] end + +@testset "New style array variable conversion and interception" begin + # Parameters, variables, and derivatives + n_comp = 2 + @parameters t, x, p[1:n_comp], q[1:n_comp] + @variables u(..)[1:n_comp] + Dt = Differential(t) + Dx = Differential(x) + Dxx = Differential(x)^2 + params = Symbolics.scalarize(reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]])) + # 1D PDE and boundary conditions + + eqs = [Dt(u(t, x)[i]) ~ p[i] * Dxx(u(t, x)[i]) for i in 1:n_comp] + + bcs = [[u(0, x)[i] ~ q[i] * cos(x), + u(t, 0)[i] ~ sin(t), + u(t, 1)[i] ~ exp(-t) * cos(1), + Dx(u(t,0)[i]) ~ 0.0] for i in 1:n_comp] + bcs_collected = reduce(vcat, bcs) + + # Space and time domains + domains = [t ∈ Interval(0.0, 1.0), + x ∈ Interval(0.0, 1.0)] + + # PDE system + + @named pdesys = PDESystem(eqs, bcs_collected, domains, [t, x], [u(t, x)[i] for i in 1:n_comp], Symbolics.scalarize(params)) + + + # Method of lines discretization + dx = 0.1 + order = 2 + discretization = MOLFiniteDifference([x => dx], t; approx_order = order) + + # Convert the PDE problem into an ODE problem + prob = discretize(pdesys,discretization) #error occurs here + + # Solve ODE problem + sol = solve(prob, Tsit5(), saveat=0.2) + + # Test that the system is correctly constructed + varname1 = Symbol("u_Any[1]") + varname2 = Symbol("u_Any[2]") + + + vars = @variables $varname1(..), $varname2(..) + + @test sol[u(t, x)[1]] == sol[vars[1](t, x)] + @test sol[u(t, x)[2]] == sol[vars[2](t, x)] +end diff --git a/test/runtests.jl b/test/runtests.jl index ef2a80c00..5916d5acf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,12 +7,17 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin - if GROUP == "All" || GROUP == "Sol_Interface" @time @safetestset "MOLFiniteDifference Interface: Solution interface" begin include("components/solution_interface.jl") end end + if GROUP == "All" || GROUP == "MOL_Interface2" + @time @safetestset "MOLFiniteDifference Interface" begin + include("pde_systems/MOLtest2.jl") + end + end + if GROUP == "All" || GROUP == "Diffusion" @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin @@ -62,9 +67,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end end - - - if GROUP == "All" || GROUP == "Higher_Order" @time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("pde_systems/MOL_1D_HigherOrder.jl") @@ -102,12 +104,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") end end - if GROUP == "All" || GROUP == "MOL_Interface2" - @time @safetestset "MOLFiniteDifference Interface" begin - include("pde_systems/MOLtest2.jl") - end - end - if GROUP == "All" || GROUP == "2D_Diffusion" @time @safetestset "MOLFiniteDifference Interface: 2D Diffusion" begin include("pde_systems/MOL_2D_Diffusion.jl")