Symbolics.jl icon indicating copy to clipboard operation
Symbolics.jl copied to clipboard

@makearray fails to propagate shape

Open xtalax opened this issue 3 years ago • 4 comments

The following:

using Symbolics, ModelingToolkit

@parameters t, x, y
@variables u(..) udisc[1:11, 1:11](t)

Dt = Differential(t)
Dxx = Differential(t)
Dyy = Differential(t)

Igrid = CartesianIndices((1:11, 1:11))
interior = Igrid[2:end-1, 2:end-1]

α = 1.1

eq = Dt(u(t,x,y)) ~ α*(Dxx(u(t, x, y)) + Dyy(u(t, x, y)))

secondderiv(II, I1) = udisc[II-I1] - 2udisc[II] + udisc[II+I1]

function generate_Dx(Igrid, interior, udisc)
    Ix = CartesianIndex(1, 0)
    @makearray out[1:11, 1:11] begin
        out[interior] => @arrayop (i, j) secondderiv(CartesianIndex(i, j), Ix) i in 2:10 j in 2:10
    end
    return out
end

rules = [
    Dxx(u(t, x, y)) => generate_Dx(Igrid, interior, udisc), #ERROR
]

rhs = substitute(eq.rhs, rules)

Produces this Error:

ERROR: DimensionMismatch("setview did not work while assigning indices (CartesianIndex{2}[CartesianIndex(2, 2) CartesianIndex(2, 3) CartesianIndex(2, 4) CartesianIndex(2, 5) CartesianIndex(2, 6) CartesianIndex(2, 7) CartesianIndex(2, 8) CartesianIndex(2, 9) CartesianIndex(2, 10); CartesianIndex(3, 2) CartesianIndex(3, 3) CartesianIndex(3, 4) CartesianIndex(3, 5) CartesianIndex(3, 6) CartesianIndex(3, 7) CartesianIndex(3, 8) CartesianIndex(3, 9) CartesianIndex(3, 10); CartesianIndex(4, 2) CartesianIndex(4, 3) CartesianIndex(4, 4) CartesianIndex(4, 5) CartesianIndex(4, 6) CartesianIndex(4, 7) CartesianIndex(4, 8) CartesianIndex(4, 9) CartesianIndex(4, 10); CartesianIndex(5, 2) CartesianIndex(5, 3) CartesianIndex(5, 4) CartesianIndex(5, 5) CartesianIndex(5, 6) CartesianIndex(5, 7) CartesianIndex(5, 8) CartesianIndex(5, 9) CartesianIndex(5, 10); CartesianIndex(6, 2) CartesianIndex(6, 3) CartesianIndex(6, 4) CartesianIndex(6, 5) CartesianIndex(6, 6) CartesianIndex(6, 7) CartesianIndex(6, 8) CartesianIndex(6, 9) CartesianIndex(6, 10); CartesianIndex(7, 2) CartesianIndex(7, 3) CartesianIndex(7, 4) CartesianIndex(7, 5) CartesianIndex(7, 6) CartesianIndex(7, 7) CartesianIndex(7, 8) CartesianIndex(7, 9) CartesianIndex(7, 10); CartesianIndex(8, 2) CartesianIndex(8, 3) CartesianIndex(8, 4) CartesianIndex(8, 5) CartesianIndex(8, 6) CartesianIndex(8, 7) CartesianIndex(8, 8) CartesianIndex(8, 9) CartesianIndex(8, 10); CartesianIndex(9, 2) CartesianIndex(9, 3) CartesianIndex(9, 4) CartesianIndex(9, 5) CartesianIndex(9, 6) CartesianIndex(9, 7) CartesianIndex(9, 8) CartesianIndex(9, 9) CartesianIndex(9, 10); CartesianIndex(10, 2) CartesianIndex(10, 3) CartesianIndex(10, 4) CartesianIndex(10, 5) CartesianIndex(10, 6) CartesianIndex(10, 7) CartesianIndex(10, 8) CartesianIndex(10, 9) CartesianIndex(10, 10)],) to @arrayop(_[i,j] := secondderiv(CartesianIndex(i, j), CartesianIndex(1, 0))) j in 2:10, i in 2:10. LHS has size (81,) and RHS has size (9, 9) -- they need to be broadcastable.")
Stacktrace:
 [1] (::Symbolics.var"#check_assignment#133")(vw::Tuple{Matrix{CartesianIndex{2}}}, op::Symbolics.ArrayOp{Matrix{Any}})
   @ Symbolics ~/.julia/dev/Symbolics.jl/src/arrays.jl:801
 [2] (::Symbolics.var"#131#135"{Symbolics.var"#check_assignment#133"})(am::ArrayMaker{Real, Matrix{Real}}, vw::Tuple{Matrix{CartesianIndex{2}}}, op::Symbolics.ArrayOp{Matrix{Any}})
   @ Symbolics ~/.julia/dev/Symbolics.jl/src/arrays.jl:814
 [3] macro expansion
   @ ~/.julia/dev/Symbolics.jl/src/arrays.jl:145 [inlined]
 [4] generate_Dx(Igrid::CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, interior::Matrix{CartesianIndex{2}}, udisc::Symbolics.Arr{Num, 2})
   @ Main ~/Code/Julia/stencil.jl:22
 [5] top-level scope
   @ ~/Code/Julia/stencil.jl:28

caused by: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 81 and 9")
Stacktrace:
 [1] _bcs1
   @ ./broadcast.jl:498 [inlined]
 [2] _bcs
   @ ./broadcast.jl:495 [inlined]
 [3] broadcast_shape(::Tuple{Int64}, ::Tuple{Int64, Int64})
   @ Base.Broadcast ./broadcast.jl:489
 [4] (::Symbolics.var"#check_assignment#133")(vw::Tuple{Matrix{CartesianIndex{2}}}, op::Symbolics.ArrayOp{Matrix{Any}})
   @ Symbolics ~/.julia/dev/Symbolics.jl/src/arrays.jl:798
 [5] (::Symbolics.var"#131#135"{Symbolics.var"#check_assignment#133"})(am::ArrayMaker{Real, Matrix{Real}}, vw::Tuple{Matrix{CartesianIndex{2}}}, op::Symbolics.ArrayOp{Matrix{Any}})
   @ Symbolics ~/.julia/dev/Symbolics.jl/src/arrays.jl:814
 [6] macro expansion
   @ ~/.julia/dev/Symbolics.jl/src/arrays.jl:145 [inlined]
 [7] generate_Dx(Igrid::CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}, interior::Matrix{CartesianIndex{2}}, udisc::Symbolics.Arr{Num, 2})
   @ Main ~/Code/Julia/stencil.jl:22
 [8] top-level scope
   @ ~/Code/Julia/stencil.jl:28

This error appears to be caused by @arrayop incorrectly propagating the shape of interior.

xtalax avatar Jul 12 '22 17:07 xtalax

what is the shape of interior?

shashi avatar Jul 12 '22 18:07 shashi

(9, 9)

xtalax avatar Jul 12 '22 19:07 xtalax

I can't reproduce this problem

shashi avatar Aug 01 '22 17:08 shashi

Getting the same error on Symbolics v4.10 and MTK v8.18

xtalax avatar Aug 04 '22 16:08 xtalax