Invalidate system cache

MWE from https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/24

```julia
using OrdinaryDiffEq
using PreallocationTools
using LinearAlgebra
using ModelingToolkit
using SparseArrays
using BenchmarkTools

# Generate the constants
const N = 100
dx = 0.5 / N
q = dualcache(zeros(N, 1))
rR = dualcache(zeros(N, 1))
ϱu = 0
dp = 100e2
Tin = 200

# generate index lists to find the solution variables easily in solution vector
const indices_Tgas = 1:N
const indices_Tsol = indices_Tgas .+ N
const indices_mflux = indices_Tsol[end] + 1
const indices_RConcentrationSolid = indices_mflux[end] .+ (1:N)

# indices = (indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid)

p_noDualCache = (dx, # control volume length
    dp, # boundary condition: pressure differential across domain
    Tin, # boundary condition: Gas inlet temperature
    q, # caching variable for q
    rR, # caching variable for rR
    ϱu)

function DiscretizedSystemMassMatrixForm_noDualcache!(du, u, p, t)
    # N,dx, dp, Tin, q, rR, ϱu, indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid = p
    dx, dp, Tin, q, rR, ϱu = p
    Tgas = @view u[indices_Tgas]
    Tsol = @view u[indices_Tsol]
    dTgas = @view du[indices_Tgas]
    dTsol = @view du[indices_Tsol]
    mflux = @view u[indices_mflux]
    dmflux = @view du[indices_mflux]
    RConcsol = @view u[indices_RConcentrationSolid]
    dRConcsol = @view du[indices_RConcentrationSolid]

    q = 1e4 * (Tgas - Tsol)
    ϱu = mflux[1]

    rR = 5.5e6 .* exp.(.-5680 ./ (273 .+ Tsol)) .* (1 .- exp.(.-RConcsol ./ 50))

    # gas energy balance
    @inbounds begin
        dTgas[1] = ϱu * 1000 * (Tin - Tgas[1]) - q[1] * dx
    end
    @inbounds for k in 2:N
        dTgas[k] = ϱu * 1000 * (Tgas[k-1] - Tgas[k]) - q[k] * dx
    end


    # solids energy balance
    @inbounds for k in 1:N
        dTsol[k] = q[k] / (1000 * 2000)
    end

    # gas momentum balance
    @inbounds begin
        dmflux[1] = -mflux[1] + 5e-3 * sqrt(dp)
    end

    # reactant concentration equations
    @inbounds for k in 1:N
        dRConcsol[k] = -rR[k]
    end


end



# build initial conditions
u0 = zeros(N * 3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000

# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid];
Jj = [indices_Tsol; indices_RConcentrationSolid];
V = ones(length(Ii));
Msparse = sparse(Ii, Jj, V, 3 * N + 1, 3 * N + 1)


f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix = Msparse)
probDAE = ODEProblem(f_DAE, u0, (0.0, 1500), p_noDualCache)
de = modelingtoolkitize(probDAE)

# MTKized system back to numerical problem:
probDAE_MTK1 = ODEProblem(de, [], (0.0, 1500), p_noDualCache, jac = true, sparse = true)

# solve the numerical problem - works
sol_MTK1 = solve(probDAE_MTK1, QBDF())

# pretty much as fast as solving non-MTKized system with Jacobian sparsity, as in above posts.
@benchmark solve(probDAE_MTK1, QBDF())

# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)

# system de2 only has 200 states after structural simplification
probDAE_MTK2 = ODEProblem(de2, [], (0.0, 1500), p_noDualCache, jac = true, sparse = true)

# Error thrown: UndefVarError: x_{201} not defined 
sol_MTK2 = solve(probDAE_MTK2, QBDF())

# Get a new system and simplify that
de = modelingtoolkitize(probDAE)
de2 = structural_simplify(de)

# system de2 only has 200 states after structural simplification
probDAE_MTK2 = ODEProblem(de2, [], (0.0, 1500), p_noDualCache, jac = true, sparse = true)

# No error!
sol_MTK2 = solve(probDAE_MTK2, QBDF())
```

Author: ChrisRackauckas

src/systems/abstractsystem.jl

@@ -929,6 +929,7 @@ function structural_simplify(sys::AbstractSystem; simplify=false)
     sys = tearing_reassemble(state, tearing(state), simplify=simplify)
     fullstates = [map(eq->eq.lhs, observed(sys)); states(sys)]
     @set! sys.observed = topsort_equations(observed(sys), fullstates)
+    invalidate_cache!(sys)
     return sys
 end
 

test/structural_transformation/tearing.jl

@@ -212,13 +212,19 @@ ms_eqs = []
 @named ms_model = compose(_ms_model,
                           [mass])
 
+calculate_jacobian(ms_model)
+calculate_tgrad(ms_model)
+calculate_control_jacobian(ms_model)
 # Mass starts with velocity = 1
 u0 = [
       mass.s => 0.0
       mass.v => 1.0
      ]
 
 sys = structural_simplify(ms_model)
+@test sys.jac[] === ModelingToolkit.EMPTY_JAC
+@test sys.tgrad[] === ModelingToolkit.EMPTY_TGRAD
+@test sys.ctrl_jac[] === ModelingToolkit.EMPTY_JAC
 prob_complex = ODAEProblem(sys, u0, (0, 1.0))
 sol = solve(prob_complex, Tsit5())
 @test all(sol[mass.v] .== 1)