deal with getindex led dvs

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]

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)

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),

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),

test/pde_systems/MOLtest2.jl

@@ -38,72 +38,83 @@ 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
         D = 1.0 / (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 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

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")