From 77169d9da3d7392c7f45767f8b09bc95a4222cfc Mon Sep 17 00:00:00 2001 From: Harsh Singh Date: Tue, 10 Mar 2026 00:08:44 +0530 Subject: [PATCH] Forward M1/M2, DIMOFIND*, and mass matrix through radau5/radau SciML wrapper --- src/algorithms.jl | 52 +++++++++++++++++++++++++++++++++++++++++------ src/solve.jl | 18 ++++++++++++++-- 2 files changed, 62 insertions(+), 8 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index e8831a6..5484fe1 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -80,7 +80,10 @@ end SciMLBase.alg_order(alg::seulex) = 12 """ - radau(; jac_lower = nothing, jac_upper = nothing) + radau(; jac_lower = nothing, jac_upper = nothing, + M1 = nothing, M2 = nothing, + DIMOFIND1VAR = nothing, DIMOFIND2VAR = nothing, DIMOFIND3VAR = nothing, + massmatrix = nothing) Implicit Runge-Kutta (Radau IIA) method of variable order between 5 and 13. @@ -91,6 +94,12 @@ efficiently. ## Arguments - `jac_lower`: Lower bandwidth of the Jacobian for banded matrices - `jac_upper`: Upper bandwidth of the Jacobian for banded matrices +- `M1`: Number of equations with trivial structure `y_i' = y_{i+M2}` (block structure) +- `M2`: Block size for the M1 trivial equations +- `DIMOFIND1VAR`: Number of index-1 variables (for DAE index information) +- `DIMOFIND2VAR`: Number of index-2 variables (for DAE index information) +- `DIMOFIND3VAR`: Number of index-3 variables (for DAE index information) +- `massmatrix`: Mass matrix for the non-trivial block (overrides `prob.f.mass_matrix`) ## Solver Properties - Order: 5 to 13 (variable) @@ -98,14 +107,23 @@ efficiently. - Type: Implicit Runge-Kutta (Radau IIA) - Suitable for: Stiff problems and DAEs """ -struct radau{T} <: ODEInterfaceImplicitAlgorithm +struct radau{T, MM} <: ODEInterfaceImplicitAlgorithm jac_lower::T jac_upper::T + M1::Union{Nothing, Int} + M2::Union{Nothing, Int} + DIMOFIND1VAR::Union{Nothing, Int} + DIMOFIND2VAR::Union{Nothing, Int} + DIMOFIND3VAR::Union{Nothing, Int} + massmatrix::MM end SciMLBase.alg_order(alg::radau) = 13 """ - radau5(; jac_lower = nothing, jac_upper = nothing) + radau5(; jac_lower = nothing, jac_upper = nothing, + M1 = nothing, M2 = nothing, + DIMOFIND1VAR = nothing, DIMOFIND2VAR = nothing, DIMOFIND3VAR = nothing, + massmatrix = nothing) Implicit Runge-Kutta method (Radau IIA) of order 5. @@ -115,6 +133,12 @@ Radau method, often preferred when a specific order is desired or for comparison ## Arguments - `jac_lower`: Lower bandwidth of the Jacobian for banded matrices - `jac_upper`: Upper bandwidth of the Jacobian for banded matrices +- `M1`: Number of equations with trivial structure `y_i' = y_{i+M2}` (block structure) +- `M2`: Block size for the M1 trivial equations +- `DIMOFIND1VAR`: Number of index-1 variables (for DAE index information) +- `DIMOFIND2VAR`: Number of index-2 variables (for DAE index information) +- `DIMOFIND3VAR`: Number of index-3 variables (for DAE index information) +- `massmatrix`: Mass matrix for the non-trivial block (overrides `prob.f.mass_matrix`) ## Solver Properties - Order: 5 @@ -122,9 +146,15 @@ Radau method, often preferred when a specific order is desired or for comparison - Type: Implicit Runge-Kutta (Radau IIA) - Suitable for: Stiff problems and DAEs """ -struct radau5{T} <: ODEInterfaceImplicitAlgorithm +struct radau5{T, MM} <: ODEInterfaceImplicitAlgorithm jac_lower::T jac_upper::T + M1::Union{Nothing, Int} + M2::Union{Nothing, Int} + DIMOFIND1VAR::Union{Nothing, Int} + DIMOFIND2VAR::Union{Nothing, Int} + DIMOFIND3VAR::Union{Nothing, Int} + massmatrix::MM end SciMLBase.alg_order(alg::radau5) = 5 @@ -196,7 +226,17 @@ end SciMLBase.alg_order(alg::ddebdf) = 5 seulex(; jac_lower = nothing, jac_upper = nothing) = seulex(jac_lower, jac_upper) -radau(; jac_lower = nothing, jac_upper = nothing) = radau(jac_lower, jac_upper) -radau5(; jac_lower = nothing, jac_upper = nothing) = radau5(jac_lower, jac_upper) +function radau(; jac_lower = nothing, jac_upper = nothing, + M1 = nothing, M2 = nothing, + DIMOFIND1VAR = nothing, DIMOFIND2VAR = nothing, DIMOFIND3VAR = nothing, + massmatrix = nothing) + radau(jac_lower, jac_upper, M1, M2, DIMOFIND1VAR, DIMOFIND2VAR, DIMOFIND3VAR, massmatrix) +end +function radau5(; jac_lower = nothing, jac_upper = nothing, + M1 = nothing, M2 = nothing, + DIMOFIND1VAR = nothing, DIMOFIND2VAR = nothing, DIMOFIND3VAR = nothing, + massmatrix = nothing) + radau5(jac_lower, jac_upper, M1, M2, DIMOFIND1VAR, DIMOFIND2VAR, DIMOFIND3VAR, massmatrix) +end rodas(; jac_lower = nothing, jac_upper = nothing) = rodas(jac_lower, jac_upper) ddebdf(; jac_lower = nothing, jac_upper = nothing) = ddebdf(jac_lower, jac_upper) diff --git a/src/solve.jl b/src/solve.jl index addd940..a332e8a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -126,7 +126,21 @@ function DiffEqBase.__solve( ODEINTERFACE_ALIASES, ODEINTERFACE_ALIASES_REVERSED ) - if prob.f.mass_matrix != I + + # Forward DAE structure options from radau/radau5 algorithm fields + if alg isa Union{radau, radau5} + for field in (:M1, :M2, :DIMOFIND1VAR, :DIMOFIND2VAR, :DIMOFIND3VAR) + val = getfield(alg, field) + if val !== nothing + dict[field] = val + end + end + if alg.massmatrix !== nothing + dict[:MASSMATRIX] = alg.massmatrix + end + end + + if !haskey(dict, :MASSMATRIX) && prob.f.mass_matrix != I if prob.f.mass_matrix isa Matrix && isstiff dict[:MASSMATRIX] = prob.f.mass_matrix elseif !isstiff @@ -331,7 +345,7 @@ const ODEINTERFACE_OPTION_LIST = Set( :JACOBIMATRIX, :JACOBIBANDSSTRUCT, :WORKFORRHS, :WORKFORJAC, :WORKFORDEC, :WORKFORSOL, :MAXNEWTONITER, :NEWTONSTARTZERO, :NEWTONSTOPCRIT, - :DIMFIND1VAR, + :DIMOFIND1VAR, :DIMOFIND2VAR, :DIMOFIND3VAR, :MAXSTAGES, :MINSTAGES, :INITSTAGES, :STEPSIZESTRATEGY, :FREEZESSLEFT, :FREEZESSRIGHT, :ORDERDECFACTOR,