Looking for DiffEQArrayOperator replacement

General Usage
, ,
1

Hi,

I got asked to maintain a package that someone else in my research group developed, so I apologize if this question is off topic. The package I’m working on uses SciMLBase v1, but I’m looking to update to v2.

Previously there was a struct DiffEQArrayOperator (see link 1) that the guy who developed the code I’m working on used to instantiate a mass matrix and pass that to an ODEFunction. The upgrade to v2 appears to have eliminated this operator. I assume that the replacement should be a MatrixOperator (I saw this post: link 2), but I’m not really sure. I tried the MatrixOperator, but it caused a previously stable DE to go unstable, so I must have implemented something wrong. Any pointers in the correct direction would be helpful.

Spaces are added after https and the site name, because I haven’t been able to post with links included.
Link 1: https:// github. com/SciML/SciMLBase.jl/blob/eb9299515d9c8f50d08dc6406f2a929eb23d8a06/src/operators/basic_operators.jl#L91
Link 2: https:// discourse .julialang. org/t/diffeqarrayoperator-is-missing/127500/2

Here is the snippet of code that’s breaking (it can be found here: https:// github .com/byuflowlab/GXBeam.jl/blob/d9543e9130f1fd9a5be82c8ccc0b9d58c4610097/src/interfaces/diffeq.jl#L257):

else

        # system indices
        indices = SystemIndices(assembly.start, assembly.stop, static=false, expanded=false)

        # default keyword arguments
        constants = (; assembly, indices, two_dimensional, structural_damping, force_scaling,
            xpfunc, pfunc, prescribed_conditions, distributed_loads, point_masses,
            linear_velocity, angular_velocity, gravity, t=0.0)

        # set state rate vector
        du = zeros(indices.nstates) # state rate vector must be zero
        u = rand(indices.nstates) # state vector can be anything
        t = rand() # time can be anything

        # mass matrix
        TF = eltype(system)
        nx = indices.nstates
        M = zeros(TF, nx, nx)
        update_mass_matrix! = (jacob, x, p, t) -> begin
            # zero out all mass matrix entries
            M .= 0.0
            # compute mass matrix
            mass_matrix!(jacob, x, p, (; constants..., t))
            # change sign of mass matrix
            M .*= -1
            # return result
            return M
        end
        # mass_matrix = SciMLBase.DiffEqArrayOperator(M, update_func = update_mass_matrix!)
        mass_matrix = SciMLBase.MatrixOperator(M, update_func = update_mass_matrix!) #Todo: This made things unstable.... 

        # residual function
        f = (resid, u, p, t) -> dynamic_residual!(resid, du, u, p, (; constants..., t))

        # jacobian function
        update_jacobian! = (jacob, u, p, t) -> dynamic_jacobian!(jacob, du, u, p, (; constants..., t))

        # jacobian prototype
        if sparse
            jac_prototype = spzeros(TF, nx, nx)
            update_jacobian!(jac_prototype, u, p, t)
        else
            jac_prototype = nothing
        end

    end

    if isnothing(xpfunc)
        odefunc = SciMLBase.ODEFunction{true,true}(f; mass_matrix = mass_matrix,
            jac = update_jacobian!, jac_prototype = jac_prototype)
    else
        odefunc = SciMLBase.ODEFunction{true,true}(f; mass_matrix = mass_matrix)
    end
2

This piece of code is pretty suspect. Did you mean:

        update_mass_matrix! = (jacob, x, p, t) -> begin
            # zero out all mass matrix entries
            jacob .= 0.0
            # compute mass matrix
            mass_matrix!(jacob, x, p, (; constants..., t))
            # change sign of mass matrix
            jacob .*= -1
            # return result
            return jacob
        end

The previous version was relying on M and jacob having the same pointer, but that won’t necessarily always be the case.