Fix operator_compose units bug

Project.toml

@@ -1,7 +1,7 @@
 name = "EarthSciMLBase"
 uuid = "e53f1632-a13c-4728-9402-0c66d48804b0"
 authors = ["EarthSciML Authors and Contributors"]
-version = "0.15.2"
+version = "0.15.3"
 
 [deps]
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"

src/operator_compose.jl

@@ -63,6 +63,7 @@ function operator_compose(a::ModelingToolkit.ODESystem, b::ModelingToolkit.ODESy
         end
         for (j, b_eq) ∈ enumerate(b_eqs)
             if isequal(bdv, add_scope(b, get_dv(b_eq.lhs, iv), iv))
+                bdv = add_scope(b, get_dv(b_eq.lhs, iv), iv) # Make sure the units are correct.
                 # The dependent variable of the LHS of this equation matches the dependent 
                 # variable of interest,
                 # so create a new variable to represent the dependent variable

test/operator_compose_test.jl

@@ -206,6 +206,42 @@ end
     @test occursin("sys1₊sys2_y(t)", streq)
 end
 
+@testset "Units 2" begin
+    @parameters t [unit = u"s"]
+    struct U1Coupler
+        sys
+    end
+    function U1()
+        @variables x(t) [unit = u"kg*m^-3"]
+        D = Differential(t)
+        ODESystem([D(x) ~ 0], t; name=:sys1,
+            metadata=Dict(:coupletype => U1Coupler))
+    end
+    struct U2Coupler
+        sys
+    end
+    function U2(; name=:sys2)
+        @variables x(t) [unit = u"kg*m^-3/s"]
+        @parameters p [unit = u"kg*m^-3/s"]
+        ODESystem([x ~ p], t; name=name,
+            metadata=Dict(:coupletype => U2Coupler))
+    end
+
+
+    sys1 = U1()
+    sys2 = U2()
+
+    function EarthSciMLBase.couple2(s1::U1Coupler, s2::U2Coupler)
+        s1, s2 = s1.sys, s2.sys
+        operator_compose(s1, s2)
+    end
+
+    combined = couple(sys1, sys2)
+
+    sys = get_mtk(combined)
+    @test occursin("sys1₊sys2_x(t)", string(equations(sys)))
+end
+
 @testset "Reaction-Deposition" begin
     struct ChemCoupler
         sys