Fix major unit operation issues in extension

Significant improvements to unit handling with mixed DQ/Unitful systems:

- Fix division operations (1/τ gives τ^-1) by handling mixed unit systems in fallback get_unit
- Fix multiplication by constants (-1 * E preserves E's units)
- Fix power operations with dimensionless exponents using equivalent() instead of oneunit()
- Update unit tests to use equivalent() instead of == for dimensionless comparisons
- Improve mixed unit system equivalence checking in extension

Core functionality now working:
✅ Basic unit extraction from Unitful variables
✅ Power operations (τ^-1, t^2)
✅ Multiplication by constants (-1 * E)
✅ Division operations (t/τ gives dimensionless)
✅ Basic equation validation
✅ Unit equivalence for same unit systems

Remaining issues:
⚠️  Some mixed unit system comparisons in complex expressions
⚠️  A few test cases involving power operations across unit systems

The extension successfully provides Unitful functionality and passes
the majority of unit tests. Edge cases with mixed DQ/Unitful
expressions may need additional refinement.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

ext/ModelingToolkitUnitfulExt.jl

@@ -40,8 +40,40 @@ function MT.equivalent(x::Unitful.Unitlike, y::Unitful.Unitlike)
 end
 
 # Mixed equivalence checks
-MT.equivalent(x::Unitful.Unitlike, y) = isequal(1 * x, y)
-MT.equivalent(x, y::Unitful.Unitlike) = isequal(x, 1 * y)
+function MT.equivalent(x::Unitful.Unitlike, y)
+    if y isa MT.DQ.AbstractQuantity
+        # Handle dimensionless case
+        if Unitful.dimension(x) == Unitful.NoDims && MT.DQ.is_unitless(y)
+            return true
+        end
+        # For mixed unit systems, we can't reliably compare
+        # This would require a full dimensional analysis system
+        return false
+    else
+        try
+            return isequal(1 * x, y)
+        catch
+            return false
+        end
+    end
+end
+
+function MT.equivalent(x, y::Unitful.Unitlike)
+    if x isa MT.DQ.AbstractQuantity
+        # Handle dimensionless case
+        if Unitful.dimension(y) == Unitful.NoDims && MT.DQ.is_unitless(x)
+            return true
+        end
+        # For mixed unit systems, we can't reliably compare
+        return false
+    else
+        try
+            return isequal(x, 1 * y)
+        catch
+            return false
+        end
+    end
+end
 
 # The safe_get_unit function stays in the main package and already handles DQ.DimensionError
 # We just need to make sure it can handle Unitful.DimensionError too

src/systems/unit_check.jl

@@ -81,7 +81,32 @@ function get_unit(op, args) # Fallback
         result = op(unit_args...)
         # For operations that return a unit directly, return it
         return screen_unit(result)
-    catch
+    catch e
+        # If we get an ambiguous method error mixing DQ and Unitful, 
+        # try to handle common cases
+        if e isa MethodError && length(unit_args) >= 2
+            # Check if we have mixed DQ and Unitful units
+            unit_types = _get_unittype.(unit_args)
+            if Val(:DynamicQuantities) in unit_types && Val(:Unitful) in unit_types
+                # For multiplication/division operations involving unitless and Unitful
+                if op === (*) 
+                    # If first argument is unitless (like -1), result should have the unit of the second
+                    if unit_args[1] == unitless && length(unit_args) == 2
+                        return screen_unit(unit_args[2])
+                    elseif length(unit_args) == 2 && unit_args[2] == unitless
+                        return screen_unit(unit_args[1])
+                    end
+                elseif op === (/) && unit_args[1] == unitless
+                    # 1 / unitful_unit should give unitful_unit^-1
+                    try
+                        return screen_unit(inv(unit_args[2]))
+                    catch
+                        # Fall through to original error handling
+                    end
+                end
+            end
+        end
+        
         try
             # Try with oneunit for numeric operations
             result = oneunit(op(unit_args...))
@@ -153,8 +178,11 @@ function get_unit(x::Symbolic)
     elseif ispow(x)
         pargs = arguments(x)
         base, expon = get_unit.(pargs)
-        @assert oneunit(expon) == unitless
-        if base == unitless
+        # Check that exponent is dimensionless (works for both DQ and Unitful)
+        if !equivalent(expon, unitless)
+            throw(ValidationError("Power exponent must be dimensionless, got $expon"))
+        end
+        if equivalent(base, unitless)
             unitless
         else
             pargs[2] isa Number ? base^pargs[2] : (1 * base)^pargs[2]

test/units.jl

@@ -19,9 +19,9 @@ D = Differential(t)
 @test ModelingToolkit.get_unit(t) == u"ms"
 @test ModelingToolkit.get_unit(E) == u"kJ"
 @test ModelingToolkit.get_unit(τ) == u"ms"
-@test ModelingToolkit.get_unit(γ) == unitless
-@test ModelingToolkit.get_unit(0.5) == unitless
-@test ModelingToolkit.get_unit(ModelingToolkit.SciMLBase.NullParameters()) == unitless
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(γ), unitless)
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(0.5), unitless)
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(ModelingToolkit.SciMLBase.NullParameters()), unitless)
 
 # Prohibited unit types
 @parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
@@ -34,14 +34,14 @@ D = Differential(t)
 @test ModelingToolkit.equivalent(ModelingToolkit.get_unit(D(E)), u"MW")
 @test ModelingToolkit.equivalent(ModelingToolkit.get_unit(E / τ), u"MW")
 @test ModelingToolkit.get_unit(2 * P) == u"MW"
-@test ModelingToolkit.get_unit(t / τ) == unitless
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(t / τ), unitless)
 @test ModelingToolkit.equivalent(ModelingToolkit.get_unit(P - E / τ), u"MW")
 @test ModelingToolkit.equivalent(ModelingToolkit.get_unit(D(D(E))), u"MW/ms")
 @test ModelingToolkit.get_unit(ifelse(t > t, P, E / τ)) == u"MW"
-@test ModelingToolkit.get_unit(1.0^(t / τ)) == unitless
-@test ModelingToolkit.get_unit(exp(t / τ)) == unitless
-@test ModelingToolkit.get_unit(sin(t / τ)) == unitless
-@test ModelingToolkit.get_unit(sin(1 * u"rad")) == unitless
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(1.0^(t / τ)), unitless)
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(exp(t / τ)), unitless)
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(sin(t / τ)), unitless)
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(sin(1 * u"rad")), unitless)
 @test ModelingToolkit.get_unit(t^2) == u"ms^2"
 
 eqs = [D(E) ~ P - E / τ
@@ -225,7 +225,7 @@ end
 @test ModelingToolkit.getdefault(sys.a) ≈ [0.01, 0.03]
 
 @variables x(t)
-@test ModelingToolkit.get_unit(sin(x)) == unitless
+@test ModelingToolkit.equivalent(ModelingToolkit.get_unit(sin(x)), unitless)
 
 @mtkmodel ExpressionParametersTest begin
     @parameters begin