From db9dc084eccfddad71add797022e3a8884de57af Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 19 Jan 2022 22:56:34 -0600
Subject: [PATCH 01/22] add internal to simulation

---
 src/PowerSimulationsDynamics.jl |  1 +
 src/base/simulation.jl          | 67 ++++++++++++++-------------------
 src/base/simulation_internal.jl | 25 ++++++++++++
 3 files changed, 54 insertions(+), 39 deletions(-)
 create mode 100644 src/base/simulation_internal.jl

diff --git a/src/PowerSimulationsDynamics.jl b/src/PowerSimulationsDynamics.jl
index 50de5984e..a8e535af2 100644
--- a/src/PowerSimulationsDynamics.jl
+++ b/src/PowerSimulationsDynamics.jl
@@ -79,6 +79,7 @@ include("base/branch_wrapper.jl")
 include("base/frequency_reference.jl")
 include("base/file_system.jl")
 include("base/simulation_model.jl")
+include("base/simulation_internal.jl")
 include("base/simulation_inputs.jl")
 include("base/perturbations.jl")
 include("base/caches.jl")
diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 89ea74f26..38ba1a1f8 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -1,30 +1,26 @@
 mutable struct Simulation{T <: SimulationModel}
     status::BUILD_STATUS
-    problem::Union{Nothing, SciMLBase.DEProblem}
-    tspan::NTuple{2, Float64}
     sys::PSY.System
     perturbations::Vector{<:Perturbation}
-    x0_init::Vector{Float64}
-    initialized::Bool
-    tstops::Vector{Float64}
-    callbacks::DiffEqBase.CallbackSet
     simulation_folder::String
     inputs::Union{Nothing, SimulationInputs}
     results::Union{Nothing, SimulationResults}
+    internal::Union{Nothing, SimulationInternal}
     console_level::Base.CoreLogging.LogLevel
     file_level::Base.CoreLogging.LogLevel
     multimachine::Bool
+    user_provided_initial_conditions::Vector{Float64}
+    initialized::Bool
 end
 
 get_system(sim::Simulation) = sim.sys
 get_simulation_inputs(sim::Simulation) = sim.inputs
-get_initial_conditions(sim::Simulation) = sim.x0_init
-get_tspan(sim::Simulation) = sim.tspan
+get_initial_conditions(sim::Simulation) = sim.internal.x0_init
+get_tspan(sim::Simulation) = sim.internal.tspan
 
 function Simulation(
     ::Type{T},
     sys::PSY.System;
-    tspan,
     initial_conditions,
     initialize_simulation,
     perturbations,
@@ -36,20 +32,17 @@ function Simulation(
 
     return Simulation{T}(
         BUILD_INCOMPLETE,
-        nothing,
-        tspan,
         sys,
         perturbations,
-        initial_conditions,
-        !initialize_simulation,
-        Vector{Float64}(),
-        DiffEqBase.CallbackSet(),
         simulation_folder,
         nothing,
         nothing,
+        nothing,
         console_level,
         file_level,
         false,
+        initial_conditions,
+        !initialize_simulation,
     )
 end
 
@@ -95,9 +88,9 @@ Builds the simulation object and conducts the indexing process. The original sys
 - `perturbations::Vector{<:Perturbation}` : Vector of Perturbations for the Simulation. Default: No Perturbations
 - `initialize_simulation::Bool` : Runs the initialization routine. If false, simulation runs based on the operation point stored in System
 - `initial_conditions::Vector{Float64}` : Allows the user to pass a vector with the initial condition values desired in the simulation. If initialize_simulation = true, these values are used as a first guess and overwritten.
-- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available: 
+- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available:
     - `ConstantFrequency` assumes that the network frequency is 1.0 per unit at all times.
-    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used. 
+    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used.
 - `system_to_file::Bool` : Default `false`. Serializes the initialized system
 - `console_level::Logging` : Default `Logging.Warn`. Sets the level of logging output to the console. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
 - `file_level::Logging` : Default `Logging.Debug`. Sets the level of logging output to file. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
@@ -152,9 +145,9 @@ Builds the simulation object and conducts the indexing process. The initial cond
 - `perturbations::Vector{<:Perturbation}` : Vector of Perturbations for the Simulation. Default: No Perturbations
 - `initialize_simulation::Bool` : Runs the initialization routine. If false, simulation runs based on the operation point stored in System
 - `initial_conditions::Vector{Float64}` : Allows the user to pass a vector with the initial condition values desired in the simulation. If initialize_simulation = true, these values are used as a first guess and overwritten.
-- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available: 
+- `frequency_reference` : Default `ReferenceBus`. Determines which frequency model is used for the network. Currently there are two options available:
     - `ConstantFrequency` assumes that the network frequency is 1.0 per unit at all times.
-    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used. 
+    - `ReferenceBus` will use the frequency state of a Dynamic Generator (rotor speed) or Dynamic Inverter (virtual speed) connected to the Reference Bus (defined in the Power Flow data) as the network frequency. If multiple devices are connected to such bus, the device with larger base power will be used as a reference. If a Voltage Source is connected to the Reference Bus, then a `ConstantFrequency` model will be used.
 - `system_to_file::Bool` : Default `false`. Serializes the initialized system
 - `console_level::Logging` : Default `Logging.Warn`. Sets the level of logging output to the console. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
 - `file_level::Logging` : Default `Logging.Debug`. Sets the level of logging output to file. Can be set to `Logging.Error`, `Logging.Warn`, `Logging.Info` or `Logging.Debug`
@@ -174,7 +167,6 @@ function Simulation(
     sim = Simulation(
         T,
         simulation_system;
-        tspan = tspan,
         initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
         initialize_simulation = get(kwargs, :initialize_simulation, true),
         simulation_folder = simulation_folder,
@@ -331,27 +323,22 @@ function _get_diffeq_problem(
     return
 end
 
-function _get_diffeq_problem(
+function _get_diffeq_function(
     sim::Simulation,
     model::SystemModel{MassMatrixModel},
     jacobian::JacobianFunctionWrapper,
 )
     simulation_inputs = get_simulation_inputs(sim)
-    sim.problem = SciMLBase.ODEProblem(
-        SciMLBase.ODEFunction{true}(
-            model,
-            mass_matrix = get_mass_matrix(simulation_inputs),
-            jac = jacobian,
-            jac_prototype = jacobian.Jv,
-            # Necessary to avoid unnecessary calculations in Rosenbrock methods
-            tgrad = (dT, u, p, t) -> dT .= false,
-        ),
-        sim.x0_init,
-        get_tspan(sim),
-        simulation_inputs,
+    val = SciMLBase.ODEFunction{true}(
+        model,
+        mass_matrix = get_mass_matrix(simulation_inputs),
+        jac = jacobian,
+        jac_prototype = jacobian.Jv,
+        # Necessary to avoid unnecessary calculations in Rosenbrock methods
+        tgrad = (dT, u, p, t) -> dT .= false,
     )
     sim.status = BUILT
-    return
+    return val
 end
 
 function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
@@ -383,7 +370,7 @@ function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
                 get_global_vars_update_pointers(sim.inputs)[GLOBAL_VAR_SYS_FREQ_INDEX] != 0
         end
         TimerOutputs.@timeit BUILD_TIMER "Pre-initialization" begin
-            _pre_initialize_simulation!(sim)
+            x0_init = _pre_initialize_simulation!(sim)
         end
         if sim.status != BUILD_FAILED
             simulation_inputs = get_simulation_inputs(sim)
@@ -392,17 +379,19 @@ function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
                     jacobian = _get_jacobian(sim)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Make Model Function" begin
-                    model = T(simulation_inputs, get_initial_conditions(sim), SimCache)
+                    model = T(simulation_inputs, x0_init, SimCache)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Initial Condition NLsolve refinement" begin
                     refine_initial_condition!(sim, model, jacobian)
                 end
+                TimerOutputs.@timeit BUILD_TIMER "Make DiffEq Problem" begin
+                    model_function = _get_model_function(sim, model, jacobian)
+                end
+
+                sim.internal = SimulationInternal(x0_init, jacobian, model, model_function)
                 TimerOutputs.@timeit BUILD_TIMER "Build Perturbations" begin
                     _build_perturbations!(sim)
                 end
-                TimerOutputs.@timeit BUILD_TIMER "Make DiffEq Problem" begin
-                    _get_diffeq_problem(sim, model, jacobian)
-                end
                 @info "Simulations status = $(sim.status)"
             catch e
                 bt = catch_backtrace()
diff --git a/src/base/simulation_internal.jl b/src/base/simulation_internal.jl
new file mode 100644
index 000000000..1ad951b5c
--- /dev/null
+++ b/src/base/simulation_internal.jl
@@ -0,0 +1,25 @@
+mutable struct SimulationInternal{T, U, V}
+    x0_init::Vector{Float64}
+    jacobian::T
+    system_model::U
+    sciml_function::V
+    tstops::Vector{Float64}
+    callbacks::DiffEqBase.CallbackSet
+    problem::Union{Nothing, SciMLBase.DEProblem}
+end
+
+function SimulationInternal(x0_init, jacobian, system_model, model_function)
+    return SimulationInternal(
+        x0_init,
+        jacobian,
+        system_model,
+        model_function,
+        Vector{Float64}(),
+        DiffEqBase.CallbackSet(),
+        nothing,
+    )
+end
+
+get_initial_conditions(internal::SimulationInternal) = internal.x0_init
+get_jacobian(internal::SimulationInternal) = internal.jacobian
+get_system_model(internal::SimulationInternal) = internal.system_model

From 612c3cda7aa12d2fb1d8ad9f62d210983c97f1de Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Wed, 19 Jan 2022 22:56:41 -0600
Subject: [PATCH 02/22] update formatter

---
 .github/workflows/formatter/Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/formatter/Project.toml b/.github/workflows/formatter/Project.toml
index 63e3e9484..1896e169e 100644
--- a/.github/workflows/formatter/Project.toml
+++ b/.github/workflows/formatter/Project.toml
@@ -6,5 +6,5 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [compat]
-JuliaFormatter = "v0.16.2"
-julia = "^1.2"
+JuliaFormatter = "v0.21"
+julia = "^1.7"

From 3cbc4374a3643c4a5be0b178f63c6df9139315bd Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 00:59:20 -0600
Subject: [PATCH 03/22] refactor simulation methods

---
 src/base/jacobian.jl                  |   8 ++
 src/base/nlsolve_wrapper.jl           |   2 +-
 src/base/simulation.jl                | 122 ++++++++++++++------------
 src/base/simulation_initialization.jl |  15 ++--
 src/utils/print.jl                    |   3 +-
 5 files changed, 81 insertions(+), 69 deletions(-)

diff --git a/src/base/jacobian.jl b/src/base/jacobian.jl
index aad83925d..78a2aaa43 100644
--- a/src/base/jacobian.jl
+++ b/src/base/jacobian.jl
@@ -96,3 +96,11 @@ function (J::JacobianFunctionWrapper)(
     JM .= gamma * LinearAlgebra.Diagonal(ones(length(x))) .- J.Jf(JM, x)
     return JM
 end
+
+function get_jacobian(
+    ::Type{T},
+    inputs::SimulationInputs,
+    x0_init::Vector{Float64},
+) where {T <: SimulationModel}
+    return JacobianFunctionWrapper(T(inputs, x0_init, JacobianCache), x0_init)
+end
diff --git a/src/base/nlsolve_wrapper.jl b/src/base/nlsolve_wrapper.jl
index d1d7ba2ac..8a19cec9c 100644
--- a/src/base/nlsolve_wrapper.jl
+++ b/src/base/nlsolve_wrapper.jl
@@ -71,6 +71,7 @@ function _convergence_check(sys_solve::NLsolveWrapper, tol::Float64, solv::Symbo
 end
 
 function refine_initial_condition!(
+    initial_guess::Vector{Float64},
     sim::Simulation,
     model::SystemModel,
     jacobian::JacobianFunctionWrapper,
@@ -84,7 +85,6 @@ function refine_initial_condition!(
 
     @debug "Start NLSolve System Run"
     converged = false
-    initial_guess = get_initial_conditions(sim)
     inputs = get_simulation_inputs(sim)
     bus_range = get_bus_range(inputs)
     powerflow_solution = deepcopy(initial_guess[bus_range])
diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 38ba1a1f8..6f2dd8db2 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -50,22 +50,20 @@ function Simulation!(
     ::Type{T},
     system::PSY.System,
     simulation_folder::String,
-    tspan::NTuple{2, Float64},
     perturbation::Perturbation;
     kwargs...,
 ) where {T <: SimulationModel}
-    return Simulation!(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
+    return Simulation!(T, system, simulation_folder, [perturbation]; kwargs...)
 end
 
 function Simulation(
     ::Type{T},
     system::PSY.System,
     simulation_folder::String,
-    tspan::NTuple{2, Float64},
     perturbation::Perturbation;
     kwargs...,
 ) where {T <: SimulationModel}
-    return Simulation(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
+    return Simulation(T, system, simulation_folder, [perturbation]; kwargs...)
 end
 
 """
@@ -73,7 +71,6 @@ end
         ::SimulationModel
         system::PowerSystems.System
         simulation_folder::String
-        tspan::NTuple{2, Float64},
         perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
         kwargs...,
     end
@@ -100,7 +97,6 @@ function Simulation!(
     ::Type{T},
     system::PSY.System,
     simulation_folder::String,
-    tspan::NTuple{2, Float64},
     perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
     kwargs...,
 ) where {T <: SimulationModel}
@@ -109,7 +105,6 @@ function Simulation!(
     sim = Simulation(
         T,
         system;
-        tspan = tspan,
         initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
         initialize_simulation = get(kwargs, :initialize_simulation, true),
         simulation_folder = simulation_folder,
@@ -130,7 +125,6 @@ end
         ::SimulationModel
         system::PowerSystems.System
         simulation_folder::String
-        tspan::NTuple{2, Float64},
         perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
         kwargs...,
     end
@@ -157,7 +151,6 @@ function Simulation(
     ::Type{T},
     system::PSY.System,
     simulation_folder::String,
-    tspan::NTuple{2, Float64},
     perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
     kwargs...,
 ) where {T <: SimulationModel}
@@ -224,14 +217,14 @@ end
 
 function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel}
     simulation_inputs = get_simulation_inputs(sim)
-    x0_init = sim.x0_init
+    x0_init = deepcopy(sim.user_provided_initial_conditions)
     if isempty(x0_init) && sim.initialized
         @warn "Initial Conditions set to flat start"
-        sim.x0_init = _get_flat_start(simulation_inputs)
+        x0_init = _get_flat_start(simulation_inputs)
     elseif isempty(x0_init) && !sim.initialized
-        sim.x0_init = _get_flat_start(simulation_inputs)
+        x0_init = _get_flat_start(simulation_inputs)
     elseif !isempty(x0_init) && sim.initialized
-        if length(sim.x0_init) != get_variable_count(simulation_inputs)
+        if length(x0_init) != get_variable_count(simulation_inputs)
             throw(
                 IS.ConflictingInputsError(
                     "The size of the provided initial state space does not match the model's state space.",
@@ -240,20 +233,20 @@ function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel
         end
     elseif !isempty(x0_init) && !sim.initialized
         @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
-        if length(sim.x0_init) != get_variable_count(simulation_inputs)
-            sim.x0_init = _get_flat_start(simulation_inputs)
+        if length(x0_init) != get_variable_count(simulation_inputs)
+            x0_init = _get_flat_start(simulation_inputs)
         end
     else
         @assert false
     end
-    return
+    return x0_init
 end
 
 function _pre_initialize_simulation!(sim::Simulation)
-    _initialize_state_space(sim)
+    x0_init = _initialize_state_space(sim)
     if sim.initialized != true
         @info("Pre-Initializing Simulation States")
-        sim.initialized = precalculate_initial_conditions!(sim)
+        sim.initialized = precalculate_initial_conditions!(x0_init, sim)
         if !sim.initialized
             error(
                 "The simulation failed to find an adequate initial guess for the initialization. Check the intialization routine.",
@@ -265,13 +258,11 @@ function _pre_initialize_simulation!(sim::Simulation)
         )
         sim.status = SIMULATION_INITIALIZED
     end
-    return
+    return x0_init
 end
 
-function _get_jacobian(sim::Simulation{T}) where {T <: SimulationModel}
-    inputs = get_simulation_inputs(sim)
-    x0_init = get_initial_conditions(sim)
-    return JacobianFunctionWrapper(T(inputs, x0_init, JacobianCache), x0_init)
+function _get_jacobian(::Simulation{T}, inputs, x0_init) where {T <: SimulationModel}
+    return get_jacobian(T, inputs, x0_init)
 end
 
 function _build_perturbations!(sim::Simulation)
@@ -292,44 +283,32 @@ function _build_perturbations!(sim::Simulation)
         callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
         tstops[ix] = pert.time
     end
-    sim.tstops = tstops
-    sim.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
+    sim.internal.tstops .= tstops
+    sim.internal.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
     return
 end
 
-function _get_diffeq_problem(
-    sim::Simulation,
+function _get_model_function(
+    ::Simulation,
     model::SystemModel{ResidualModel},
     jacobian::JacobianFunctionWrapper,
 )
-    x0 = get_initial_conditions(sim)
-    dx0 = zeros(length(x0))
-    simulation_inputs = get_simulation_inputs(sim)
-    sim.problem = SciMLBase.DAEProblem(
-        SciMLBase.DAEFunction{true}(
-            model;
-            # Currently commented for Sundials compatibility
-            #jac = jacobian,
-            tgrad = (dT, u, p, t) -> dT .= false,
-            #jac_prototype = jacobian.Jv,
-        ),
-        dx0,
-        x0,
-        get_tspan(sim),
-        simulation_inputs,
-        differential_vars = get_DAE_vector(simulation_inputs),
+    return SciMLBase.DAEFunction{true}(
+        model;
+        # Currently commented for Sundials compatibility
+        #jac = jacobian,
+        tgrad = (dT, u, p, t) -> dT .= false,
+        jac_prototype = jacobian.Jv,
     )
-    sim.status = BUILT
-    return
 end
 
-function _get_diffeq_function(
+function _get_model_function(
     sim::Simulation,
     model::SystemModel{MassMatrixModel},
     jacobian::JacobianFunctionWrapper,
 )
     simulation_inputs = get_simulation_inputs(sim)
-    val = SciMLBase.ODEFunction{true}(
+    return SciMLBase.ODEFunction{true}(
         model,
         mass_matrix = get_mass_matrix(simulation_inputs),
         jac = jacobian,
@@ -337,8 +316,6 @@ function _get_diffeq_function(
         # Necessary to avoid unnecessary calculations in Rosenbrock methods
         tgrad = (dT, u, p, t) -> dT .= false,
     )
-    sim.status = BUILT
-    return val
 end
 
 function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
@@ -376,13 +353,13 @@ function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
             simulation_inputs = get_simulation_inputs(sim)
             try
                 TimerOutputs.@timeit BUILD_TIMER "Calculate Jacobian" begin
-                    jacobian = _get_jacobian(sim)
+                    jacobian = _get_jacobian(sim, simulation_inputs, x0_init)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Make Model Function" begin
                     model = T(simulation_inputs, x0_init, SimCache)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Initial Condition NLsolve refinement" begin
-                    refine_initial_condition!(sim, model, jacobian)
+                    refine_initial_condition!(x0_init, sim, model, jacobian)
                 end
                 TimerOutputs.@timeit BUILD_TIMER "Make DiffEq Problem" begin
                     model_function = _get_model_function(sim, model, jacobian)
@@ -392,6 +369,7 @@ function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
                 TimerOutputs.@timeit BUILD_TIMER "Build Perturbations" begin
                     _build_perturbations!(sim)
                 end
+                sim.status = BUILT
                 @info "Simulations status = $(sim.status)"
             catch e
                 bt = catch_backtrace()
@@ -450,19 +428,46 @@ function _prog_meter_enabled()
            (get(ENV, "RUNNING_PSID_TESTS", nothing) != "true")
 end
 
-function _execute!(sim::Simulation, solver; kwargs...)
+function _get_diff_eq_problem(
+    internal::SimulationInternal{T, SystemModel{ResidualModel, W}, V},
+    simulation_inputs::SimulationInputs,
+    tspan,
+) where {T, V, W}
+    x0 = get_initial_conditions(internal)
+    dx0 = similar(x0)
+    return SciMLBase.DAEProblem(
+        internal.sciml_function,
+        dx0,
+        x0,
+        tspan,
+        simulation_inputs,
+        differential_vars = get_DAE_vector(simulation_inputs),
+    )
+end
+
+function _get_diff_eq_problem(
+    internal::SimulationInternal{T, SystemModel{MassMatrixModel, W}, V},
+    simulation_inputs::SimulationInputs,
+    tspan,
+) where {T, V, W}
+    x0 = get_initial_conditions(internal)
+    return SciMLBase.ODEProblem(internal.sciml_function, x0, tspan, simulation_inputs)
+end
+
+function _execute!(sim::Simulation, solver, tspan; kwargs...)
     @debug "status before execute" sim.status
     simulation_pre_step!(sim, get(kwargs, :reset_simulation, false))
     sim.status = SIMULATION_STARTED
     time_log = Dict{Symbol, Any}()
+    problem = _get_diff_eq_problem(sim.internal, get_simulation_inputs(sim), tspan)
     solution,
     time_log[:timed_solve_time],
     time_log[:solve_bytes_alloc],
     time_log[:sec_in_gc] = @timed SciMLBase.solve(
-        sim.problem,
+        problem,
         solver;
-        callback = sim.callbacks,
-        tstops = sim.tstops,
+        callback = sim.internal.callbacks,
+        tstops = sim.internal.tstops,
         progress = get(kwargs, :enable_progress_bar, _prog_meter_enabled()),
         progress_steps = 1,
         kwargs...,
@@ -484,7 +489,8 @@ end
 """
     execute!(
         sim::Simulation,
-        solver;
+        solver,
+        tspan;
         kwargs...
     )
 
@@ -496,11 +502,11 @@ Solves the time-domain dynamic simulation model.
 - `enable_progress_bar::Bool` : Default: `true`. Enables progress bar for the integration routine.
 - Additional solver keyword arguments can be included. See [Common Solver Options](https://diffeq.sciml.ai/stable/basics/common_solver_opts/) in the `DifferentialEquations.jl` documentation for more details.
 """
-function execute!(sim::Simulation, solver; kwargs...)
+function execute!(sim::Simulation, solver, tspan; kwargs...)
     logger = configure_logging(sim, "a"; kwargs...)
     try
         Logging.with_logger(logger) do
-            _execute!(sim, solver; kwargs...)
+            _execute!(sim, solver, tspan; kwargs...)
         end
     catch e
         @error "Execution failed" exception = (e, catch_backtrace())
diff --git a/src/base/simulation_initialization.jl b/src/base/simulation_initialization.jl
index 27fbc3df1..d094dadc6 100644
--- a/src/base/simulation_initialization.jl
+++ b/src/base/simulation_initialization.jl
@@ -124,35 +124,34 @@ end
 
 # Default implementation for both models. This implementation is to future proof if there is
 # a divergence between the required build methods
-function _calculate_initial_guess!(sim::Simulation)
+function _calculate_initial_guess!(x0_init::Vector{Float64}, sim::Simulation)
     inputs = get_simulation_inputs(sim)
     @assert sim.status == BUILD_INCOMPLETE
     while sim.status == BUILD_INCOMPLETE
         @debug "Start state intialization routine"
         TimerOutputs.@timeit BUILD_TIMER "Power Flow solution" begin
-            sim.status = _power_flow_solution!(sim.x0_init, get_system(sim), inputs)
+            sim.status = _power_flow_solution!(x0_init, get_system(sim), inputs)
         end
         TimerOutputs.@timeit BUILD_TIMER "Initialize Static Injectors" begin
             sim.status = _initialize_static_injection!(inputs)
         end
         TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Injectors" begin
-            sim.status =
-                _initialize_dynamic_injection!(sim.x0_init, inputs, get_system(sim))
+            sim.status = _initialize_dynamic_injection!(x0_init, inputs, get_system(sim))
         end
         if has_dyn_lines(inputs)
             TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Branches" begin
-                sim.status = _initialize_dynamic_branches!(sim.x0_init, inputs)
+                sim.status = _initialize_dynamic_branches!(x0_init, inputs)
             end
         else
             @debug "No Dynamic Branches in the system"
         end
-        sim.status = check_valid_values(sim.x0_init, inputs)
+        sim.status = check_valid_values(x0_init, inputs)
     end
     return
 end
 
-function precalculate_initial_conditions!(sim::Simulation)
-    _calculate_initial_guess!(sim)
+function precalculate_initial_conditions!(x0_init::Vector{Float64}, sim::Simulation)
+    _calculate_initial_guess!(x0_init, sim)
     return sim.status != BUILD_FAILED
 end
 
diff --git a/src/utils/print.jl b/src/utils/print.jl
index 64bacb2df..18f6f5d36 100644
--- a/src/utils/print.jl
+++ b/src/utils/print.jl
@@ -69,8 +69,7 @@ function show_simulation_table(
         "Simulation Type" val_model
         "Initialized?" val_initialized
         "Multimachine system?" val_multimachine
-        "Time Span" string(sim.tspan)
-        "Number of States" string(length(sim.x0_init))
+        "Number of States" string(length(sim.internal.x0_init))
         "Number of Perturbations" string(length(sim.perturbations))
     ]
     PrettyTables.pretty_table(

From 380bf011e23a622b25d2cde94683b1852ce8c5f6 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 12:04:56 -0600
Subject: [PATCH 04/22] fix printing

---
 src/utils/print.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utils/print.jl b/src/utils/print.jl
index 18f6f5d36..310779caf 100644
--- a/src/utils/print.jl
+++ b/src/utils/print.jl
@@ -69,7 +69,7 @@ function show_simulation_table(
         "Simulation Type" val_model
         "Initialized?" val_initialized
         "Multimachine system?" val_multimachine
-        "Number of States" string(length(sim.internal.x0_init))
+        "Number of States" (!isnothing(sim.internal) ? string(length(sim.internal.x0_init)) : "Not built")
         "Number of Perturbations" string(length(sim.perturbations))
     ]
     PrettyTables.pretty_table(

From b6524f8ed53f7a1452edd79aa3409181c3049f57 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 12:05:02 -0600
Subject: [PATCH 05/22] fix base tests

---
 src/base/simulation_initialization.jl | 18 ++++---
 test/test_base.jl                     | 72 ++++++++++-----------------
 2 files changed, 38 insertions(+), 52 deletions(-)

diff --git a/src/base/simulation_initialization.jl b/src/base/simulation_initialization.jl
index d094dadc6..4a7216bb5 100644
--- a/src/base/simulation_initialization.jl
+++ b/src/base/simulation_initialization.jl
@@ -169,10 +169,16 @@ function read_initial_conditions(sim::Simulation)
     for bus in PSY.get_components(PSY.Bus, system)
         bus_n = PSY.get_number(bus)
         bus_ix = get_lookup(simulation_inputs)[bus_n]
-        V_R[bus_n] = sim.x0_init[bus_ix]
-        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
-        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
-        θ[bus_n] = atan(sim.x0_init[bus_ix + bus_size] / sim.x0_init[bus_ix])
+        V_R[bus_n] = get_initial_conditions(sim)[bus_ix]
+        V_I[bus_n] = get_initial_conditions(sim)[bus_ix + bus_size]
+        Vm[bus_n] = sqrt(
+            get_initial_conditions(sim)[bus_ix]^2 +
+            get_initial_conditions(sim)[bus_ix + bus_size]^2,
+        )
+        θ[bus_n] = atan(
+            get_initial_conditions(sim)[bus_ix + bus_size] /
+            get_initial_conditions(sim)[bus_ix],
+        )
     end
     results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
     for device in get_dynamic_injectors(simulation_inputs)
@@ -181,7 +187,7 @@ function read_initial_conditions(sim::Simulation)
         global_index = get_global_index(device)
         x0_device = Dict{Symbol, Float64}()
         for s in states
-            x0_device[s] = sim.x0_init[global_index[s]]
+            x0_device[s] = get_initial_conditions(sim)[global_index[s]]
         end
         results[name] = x0_device
     end
@@ -193,7 +199,7 @@ function read_initial_conditions(sim::Simulation)
             global_index = get_global_index(br)
             x0_br = Dict{Symbol, Float64}()
             for s in states
-                x0_br[s] = sim.x0_init[global_index[s]]
+                x0_br[s] = get_initial_conditions(sim)[global_index[s]]
             end
             printed_name = "Line " * name
             results[printed_name] = x0_br
diff --git a/test/test_base.jl b/test/test_base.jl
index b50c48091..0ac36c579 100644
--- a/test/test_base.jl
+++ b/test/test_base.jl
@@ -27,7 +27,6 @@ end
             ResidualModel,
             omib_sys, #system
             path1,
-            (0.0, 30.0), #time span
             Ybus_change;
             system_to_file = true,
         )
@@ -61,7 +60,6 @@ end
             ResidualModel,
             omib_sys, #system
             path2,
-            (0.0, 30.0), #time span
             Ybus_change;
             system_to_file = true,
         )
@@ -110,7 +108,7 @@ end
     end
 
     # Tests for all Dynamic Lines
-    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir(), (0.0, 10.0))
+    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir())
     @test sim.status == PSID.BUILT
     sim_inputs = sim.inputs
     DAE_vector = PSID.get_DAE_vector(sim_inputs)
@@ -136,7 +134,7 @@ end
     # Tests for dynamic lines with b = 0
     set_b!(dyn_branch12, (from = 0.0, to = 0.0))
     set_b!(dyn_branch23, (from = 0.0, to = 0.0))
-    sim = Simulation(ResidualModel, threebus_sys_dyns, pwd(), (0.0, 10.0))
+    sim = Simulation(ResidualModel, threebus_sys_dyns, pwd())
     @test sim.status == PSID.BUILT
     sim_inputs = sim.inputs
     DAE_vector = PSID.get_DAE_vector(sim_inputs)
@@ -171,22 +169,16 @@ end
         ResidualModel,
         sys,
         pwd(),
-        (0.0, 20.0),
         BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1");
         initialize_simulation = false,
         initial_conditions = x0_test,
     )
-    @test LinearAlgebra.norm(sim.x0_init - x0_test) <= 1e-6
+    @test LinearAlgebra.norm(PSID.get_initial_conditions(sim) - x0_test) <= 1e-6
     #Initialize System normally
-    sim_normal = Simulation(
-        ResidualModel,
-        sys,
-        pwd(),
-        (0.0, 20.0),
-        BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"),
-    )
+    sim_normal =
+        Simulation(ResidualModel, sys, pwd(), BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"))
     #Save states without generator at bus 2
-    x0 = sim_normal.x0_init
+    x0 = PSID.get_initial_conditions(sim_normal)
     x0_no_gen = vcat(x0[1:6], x0[13:end])
     #Make Generator 2 unavailable and transform bus into PQ bus
     gen = PSY.get_component(ThermalStandard, sys, "generator-102-1")
@@ -197,41 +189,35 @@ end
     sim_trip_gen = Simulation(
         ResidualModel,
         sys,
-        pwd(),
-        (0.0, 20.0);
+        pwd();
         initialize_simulation = false,
         initial_conditions = x0_no_gen,
     )
-    @test LinearAlgebra.norm(sim_trip_gen.x0_init - x0_no_gen) <= 1e-6
+    @test LinearAlgebra.norm(PSID.get_initial_conditions(sim_trip_gen) - x0_no_gen) <= 1e-6
     #Create Simulation without Gen 2 at steady state
-    sim_normal_no_gen = Simulation(
-        ResidualModel,
-        sys,
-        pwd(),
-        (0.0, 20.0),
-        BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"),
-    )
-    @test length(sim_normal_no_gen.x0_init) == 17
+    sim_normal_no_gen =
+        Simulation(ResidualModel, sys, pwd(), BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"))
+    @test length(PSID.get_initial_conditions(sim_normal_no_gen)) == 17
     #Ignore Initial Conditions without passing initialize_simulation = false
-    sim_ignore_init =
-        Simulation(ResidualModel, sys, pwd(), (0.0, 20.0); initial_conditions = x0_no_gen)
-    @test LinearAlgebra.norm(sim_ignore_init.x0_init - sim_normal_no_gen.x0_init) <= 1e-6
+    sim_ignore_init = Simulation(ResidualModel, sys, pwd(); initial_conditions = x0_no_gen)
+    @test LinearAlgebra.norm(
+        PSID.get_initial_conditions(sim_ignore_init) -
+        PSID.get_initial_conditions(sim_normal_no_gen),
+    ) <= 1e-6
     #Pass wrong vector size
     x0_wrong = zeros(20)
     @test_logs (:error, "Build failed") match_mode = :any Simulation(
         ResidualModel,
         sys,
-        pwd(),
-        (0.0, 20.0);
+        pwd();
         initialize_simulation = false,
         initial_conditions = x0_wrong,
     )
     #Flat start initialization
-    sim_flat =
-        Simulation(ResidualModel, sys, pwd(), (0.0, 20.0); initialize_simulation = false)
+    sim_flat = Simulation(ResidualModel, sys, pwd(); initialize_simulation = false)
     x0_flat = zeros(17)
     x0_flat[1:3] .= 1.0
-    @test LinearAlgebra.norm(sim_flat.x0_init - x0_flat) <= 1e-6
+    @test LinearAlgebra.norm(PSID.get_initial_conditions(sim_flat) - x0_flat) <= 1e-6
 end
 
 @testset "Test Network Kirchoff Calculation" begin
@@ -425,7 +411,7 @@ end
     end
 
     # Tests for all Dynamic Lines
-    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir(), (0.0, 10.0))
+    sim = Simulation(ResidualModel, threebus_sys_dyns, mktempdir())
     global_index = PSID.make_global_state_map(sim.inputs)
     @test_throws ErrorException PSID.get_state_from_ix(global_index, 40)
     @test ("V_2", :R) == PSID.get_state_from_ix(global_index, 2)
@@ -443,8 +429,7 @@ end
     sim = Simulation(
         ResidualModel,
         sys, #system
-        mktempdir(),
-        (0.0, 20.0), #time span
+        mktempdir();
         # Not initialized to speed up the test
         initialize_simulation = false,
     )
@@ -454,8 +439,7 @@ end
     sim = Simulation(
         ResidualModel,
         sys, #system
-        mktempdir(),
-        (0.0, 20.0),
+        mktempdir();
         # Not initialized to speed up the test
         initialize_simulation = false,
         frequency_reference = ConstantFrequency, #time span
@@ -473,8 +457,7 @@ end
     @test_logs (:error, "Build failed") match_mode = :any Simulation(
         ResidualModel,
         sys, #system
-        mktempdir(),
-        (0.0, 20.0),
+        mktempdir();
         # Not initialized to speed up the test
         initialize_simulation = false,
         frequency_reference = ConstantFrequency, #time span
@@ -487,8 +470,7 @@ end
     sim = Simulation(
         ResidualModel,
         omib_sys_copy, #system
-        mktempdir(),
-        (0.0, 30.0); #time span
+        mktempdir();
         all_lines_dynamic = true,
     )
     @test sim.status == PSID.BUILT
@@ -499,8 +481,7 @@ end
     sim = Simulation(
         ResidualModel,
         omib_sys_copy, #system
-        mktempdir(),
-        (0.0, 30.0); #time span
+        mktempdir();
         all_branches_dynamic = true,
     )
     @test sim.status == PSID.BUILT
@@ -511,8 +492,7 @@ end
     sim = Simulation!(
         ResidualModel,
         omib_sys_copy, #system
-        mktempdir(),
-        (0.0, 30.0); #time span
+        mktempdir(); #time span
         all_lines_dynamic = true,
     )
     @test sim.status == PSID.BUILT

From 83b84e925a4a9ecc8c8922f7cc76a40a92e8ebdd Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 12:18:42 -0600
Subject: [PATCH 06/22] use getter functions

---
 src/base/small_signal.jl                 |  7 ++++---
 src/post_processing/post_proc_results.jl | 12 ++++++------
 test/utils/get_results.jl                |  8 ++++----
 3 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/base/small_signal.jl b/src/base/small_signal.jl
index 8a1c3b257..ea1e70de2 100644
--- a/src/base/small_signal.jl
+++ b/src/base/small_signal.jl
@@ -154,13 +154,14 @@ function _get_participation_factors(
     return participation_factors
 end
 
-function small_signal_analysis(sim::Simulation; kwargs...)
+function small_signal_analysis(sim::Simulation{T}; kwargs...) where T <: SimulationModel
     #simulation_pre_step!(sim, get(kwargs, :reset_simulation, false), Real)
     #sim.status = CONVERTED_FOR_SMALL_SIGNAL
     inputs = get_simulation_inputs(sim)
     mass_matrix = get_mass_matrix(inputs)
-    x_eval = get(kwargs, :operating_point, sim.x0_init)
-    jacwrapper = _get_jacobian(sim)
+    x_eval = get(kwargs, :operating_point, PSID.get_initial_conditions(sim))
+    x0_init = get_initial_conditions(sim)
+    jacwrapper = get_jacobian(T, inputs, x0_init)
     jacwrapper(x_eval)
     jacobian = jacwrapper.Jv
     diff_states = get_DAE_vector(inputs)
diff --git a/src/post_processing/post_proc_results.jl b/src/post_processing/post_proc_results.jl
index b78c862d1..c8d83e96e 100644
--- a/src/post_processing/post_proc_results.jl
+++ b/src/post_processing/post_proc_results.jl
@@ -171,7 +171,7 @@ end
 """
     show_states_initial_value(res::SimulationResults)
 
-Function to print initial states. 
+Function to print initial states.
 
 # Arguments
 
@@ -235,7 +235,7 @@ end
 """
     show_states_initial_value(sim::Simulation)
 
-Function to print initial states. 
+Function to print initial states.
 
 # Arguments
 
@@ -256,8 +256,8 @@ function show_states_initial_value(sim::Simulation)
         println("====================")
         bus_n = PSY.get_number(bus)
         bus_ix = PSID.get_lookup(sim.inputs)[bus_n]
-        V_R = sim.x0_init[bus_ix]
-        V_I = sim.x0_init[bus_ix + bus_size]
+        V_R = PSID.get_initial_conditions(sim)[bus_ix]
+        V_I = PSID.get_initial_conditions(sim)[bus_ix + bus_size]
         Vm = sqrt(V_R^2 + V_I^2)
         θ = angle(V_R + V_I * 1im)
         print("Vm ", round(Vm, digits = 4), "\n")
@@ -273,7 +273,7 @@ function show_states_initial_value(sim::Simulation)
         println("====================")
         global_index = global_state_map[name]
         for s in states
-            print(s, " ", round(sim.x0_init[global_index[s]], digits = 4), "\n")
+            print(s, " ", round(PSID.get_initial_conditions(sim)[global_index[s]], digits = 4), "\n")
         end
         println("====================")
     end
@@ -290,7 +290,7 @@ function show_states_initial_value(sim::Simulation)
             global_index = global_state_map[name]
             x0_br = Dict{Symbol, Float64}()
             for (i, s) in enumerate(states)
-                print(s, " ", round(sim.x0_init[global_index[s]], digits = 5), "\n")
+                print(s, " ", round(PSID.get_initial_conditions(sim)[global_index[s]], digits = 5), "\n")
             end
             println("====================")
         end
diff --git a/test/utils/get_results.jl b/test/utils/get_results.jl
index b49352e82..6a79dfcff 100644
--- a/test/utils/get_results.jl
+++ b/test/utils/get_results.jl
@@ -8,8 +8,8 @@ function get_init_values_for_comparison(sim::Simulation)
     for bus in PSY.get_components(PSY.Bus, system)
         bus_n = PSY.get_number(bus)
         bus_ix = PSID.get_lookup(sim.inputs)[bus_n]
-        V_R[bus_ix] = sim.x0_init[bus_ix]
-        V_I[bus_ix] = sim.x0_init[bus_ix + bus_size]
+        V_R[bus_ix] = PSID.get_initial_conditions(sim)[bus_ix]
+        V_I[bus_ix] = PSID.get_initial_conditions(sim)[bus_ix + bus_size]
         Vm[bus_ix] = sqrt(V_R[bus_ix]^2 + V_I[bus_ix]^2)
         θ[bus_ix] = angle(V_R[bus_ix] + V_I[bus_ix] * 1im)
     end
@@ -21,7 +21,7 @@ function get_init_values_for_comparison(sim::Simulation)
         global_index = PSID.get_global_index(device)
         x0_device = Vector{Float64}(undef, length(states))
         for (i, s) in enumerate(states)
-            x0_device[i] = sim.x0_init[global_index[s]]
+            x0_device[i] = PSID.get_initial_conditions(sim)[global_index[s]]
         end
         results[name] = x0_device
     end
@@ -31,7 +31,7 @@ function get_init_values_for_comparison(sim::Simulation)
         global_index = PSID.get_global_index(br)
         x0_br = Vector{Float64}(undef, length(states))
         for (i, s) in enumerate(states)
-            x0_br[i] = sim.x0_init[global_index[s]]
+            x0_br[i] = PSID.get_initial_conditions(sim)[global_index[s]]
         end
         printed_name = "Line " * name
         results[printed_name] = x0_br

From 2eef9445a35cb4ed0828980bdc79fe907ca688d4 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 12:18:55 -0600
Subject: [PATCH 07/22] change location of tspan

---
 test/test_case01_OMIB.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/test/test_case01_OMIB.jl b/test/test_case01_OMIB.jl
index 26c478d86..63582759d 100644
--- a/test/test_case01_OMIB.jl
+++ b/test/test_case01_OMIB.jl
@@ -28,7 +28,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             omib_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change,
         )
 
@@ -51,7 +50,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test01_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -97,7 +96,6 @@ end
             MassMatrixModel,
             omib_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change,
         )
 
@@ -132,7 +130,7 @@ end
         @test LinearAlgebra.norm(eigs - test01_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 

From a018ddff975542280807e81ae77ee0d6764c6eed Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 14:40:15 -0600
Subject: [PATCH 08/22] remove DiffEqBase dependency

---
 src/PowerSimulationsDynamics.jl | 1 -
 src/base/simulation.jl          | 8 ++++----
 src/base/simulation_internal.jl | 4 ++--
 3 files changed, 6 insertions(+), 7 deletions(-)

diff --git a/src/PowerSimulationsDynamics.jl b/src/PowerSimulationsDynamics.jl
index a8e535af2..25cf7aa2e 100644
--- a/src/PowerSimulationsDynamics.jl
+++ b/src/PowerSimulationsDynamics.jl
@@ -48,7 +48,6 @@ import InfrastructureSystems
 import SciMLBase
 import DataStructures: OrderedDict
 import Random
-import DiffEqBase
 import ForwardDiff
 import SparseArrays
 import LinearAlgebra
diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 6f2dd8db2..5b1ee7813 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -269,22 +269,22 @@ function _build_perturbations!(sim::Simulation)
     @info "Attaching Perturbations"
     if isempty(sim.perturbations)
         @debug "The simulation has no perturbations"
-        return DiffEqBase.CallbackSet(), [0.0]
+        return SciMLBase.CallbackSet(), [0.0]
     end
     inputs = get_simulation_inputs(sim)
     perturbations = sim.perturbations
     perturbations_count = length(perturbations)
-    callback_vector = Vector{DiffEqBase.DiscreteCallback}(undef, perturbations_count)
+    callback_vector = Vector{SciMLBase.DiscreteCallback}(undef, perturbations_count)
     tstops = Vector{Float64}(undef, perturbations_count)
     for (ix, pert) in enumerate(perturbations)
         @debug pert
         condition = (x, t, integrator) -> t in [pert.time]
         affect = get_affect(inputs, get_system(sim), pert)
-        callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
+        callback_vector[ix] = SciMLBase.DiscreteCallback(condition, affect)
         tstops[ix] = pert.time
     end
     sim.internal.tstops .= tstops
-    sim.internal.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
+    sim.internal.callbacks = SciMLBase.CallbackSet((), tuple(callback_vector...))
     return
 end
 
diff --git a/src/base/simulation_internal.jl b/src/base/simulation_internal.jl
index 1ad951b5c..990be5de0 100644
--- a/src/base/simulation_internal.jl
+++ b/src/base/simulation_internal.jl
@@ -4,7 +4,7 @@ mutable struct SimulationInternal{T, U, V}
     system_model::U
     sciml_function::V
     tstops::Vector{Float64}
-    callbacks::DiffEqBase.CallbackSet
+    callbacks::SciMLBase.CallbackSet
     problem::Union{Nothing, SciMLBase.DEProblem}
 end
 
@@ -15,7 +15,7 @@ function SimulationInternal(x0_init, jacobian, system_model, model_function)
         system_model,
         model_function,
         Vector{Float64}(),
-        DiffEqBase.CallbackSet(),
+        SciMLBase.CallbackSet(),
         nothing,
     )
 end

From a8cedd98f5165e0272d72740965c3bf32f15b979 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 14:58:45 -0600
Subject: [PATCH 09/22] fix assignment of tstops

---
 src/base/simulation.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 5b1ee7813..0d68967f3 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -283,7 +283,7 @@ function _build_perturbations!(sim::Simulation)
         callback_vector[ix] = SciMLBase.DiscreteCallback(condition, affect)
         tstops[ix] = pert.time
     end
-    sim.internal.tstops .= tstops
+    sim.internal.tstops = tstops
     sim.internal.callbacks = SciMLBase.CallbackSet((), tuple(callback_vector...))
     return
 end

From 559b454df5087e0113522b0883d71fd4df746788 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 14:59:02 -0600
Subject: [PATCH 10/22] add more tests

---
 test/test_base.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/test_base.jl b/test/test_base.jl
index 0ac36c579..901656b22 100644
--- a/test/test_base.jl
+++ b/test/test_base.jl
@@ -38,6 +38,8 @@ end
         @test isa(dic_init_conds, Dict)
         dic_control_refs = get_setpoints(sim)
         @test isa(dic_control_refs, Dict)
+        @test !isempty(sim.internal.tstops)
+        @test !isempty(sim.internal.callbacks)
 
         o_system = System(joinpath(path1, "input_system.json"))
         for b in get_components(Bus, o_system)

From 38efa0fd965ba67a4c231e2478bae04e6f97bd9b Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 14:59:16 -0600
Subject: [PATCH 11/22] remove DiffEqBase from Project toml

---
 Project.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0d5d89ac0..02d87c74b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,6 @@ version = "0.8.5"
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
@@ -22,7 +21,6 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 DataStructures = "~0.18"
-DiffEqBase = "^6.8"
 DocStringExtensions = "~0.8.6"
 ForwardDiff = "~v0.10"
 InfrastructureSystems = "^1.14"

From 46965495af04589b696f79108fe01617e54a4ff6 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:05:23 -0600
Subject: [PATCH 12/22] update docs

---
 docs/src/execute.md           | 4 ++--
 docs/src/quick_start_guide.md | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/execute.md b/docs/src/execute.md
index 0cdf11e6b..7c21159ae 100644
--- a/docs/src/execute.md
+++ b/docs/src/execute.md
@@ -75,7 +75,7 @@ Any solver option available in `DifferentialEquations.jl` can be passed as keywo
 Most common solver options used are `dtmax` to control the maximum dt for adaptive timestepping. `abstol` and `reltol` are also commonly used to control the tolerance in the adaptive timestepping. `saveat` is also used to store the results at a specified time stamps. For example, the following code is valid to further specify your solver options:
 
 ```julia
-execute!(sim, IDA(), dtmax = 0.01, abstol = 1e-9, reltol = 1e-6, saveat = 0.01)
+execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.01, abstol = 1e-9, reltol = 1e-6, saveat = 0.01)
 ```
 
-In addition, the keyword argument `enable_progress_bar = false` can be used to disable the progress bar.
\ No newline at end of file
+In addition, the keyword argument `enable_progress_bar = false` can be used to disable the progress bar.
diff --git a/docs/src/quick_start_guide.md b/docs/src/quick_start_guide.md
index d1c92a8a4..c4cb643b4 100644
--- a/docs/src/quick_start_guide.md
+++ b/docs/src/quick_start_guide.md
@@ -64,7 +64,7 @@ This means that the state ω of the generator at bus 102, participates 50% in ei
 ## Execute the simulation
 
 ```@repl quick_start_guide
-    execute!(sim, IDA(), dtmax = 0.02, saveat = 0.02, enable_progress_bar = false)
+    execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02, saveat = 0.02, enable_progress_bar = false)
 ```
 
 ## Make a plot of the results

From 437cbd8a9187a38a2bb01e0d286e2b7c42cd760a Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:10:27 -0600
Subject: [PATCH 13/22] update small signal code

---
 src/base/small_signal.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/base/small_signal.jl b/src/base/small_signal.jl
index ea1e70de2..20bfdd2a7 100644
--- a/src/base/small_signal.jl
+++ b/src/base/small_signal.jl
@@ -154,7 +154,7 @@ function _get_participation_factors(
     return participation_factors
 end
 
-function small_signal_analysis(sim::Simulation{T}; kwargs...) where T <: SimulationModel
+function small_signal_analysis(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
     #simulation_pre_step!(sim, get(kwargs, :reset_simulation, false), Real)
     #sim.status = CONVERTED_FOR_SMALL_SIGNAL
     inputs = get_simulation_inputs(sim)

From f758ccff2668accf4636b3052b409464762a8415 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:10:37 -0600
Subject: [PATCH 14/22] remove predicate

---
 src/post_processing/post_proc_results.jl | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/post_processing/post_proc_results.jl b/src/post_processing/post_proc_results.jl
index c8d83e96e..3eb7f5adf 100644
--- a/src/post_processing/post_proc_results.jl
+++ b/src/post_processing/post_proc_results.jl
@@ -256,8 +256,8 @@ function show_states_initial_value(sim::Simulation)
         println("====================")
         bus_n = PSY.get_number(bus)
         bus_ix = PSID.get_lookup(sim.inputs)[bus_n]
-        V_R = PSID.get_initial_conditions(sim)[bus_ix]
-        V_I = PSID.get_initial_conditions(sim)[bus_ix + bus_size]
+        V_R = get_initial_conditions(sim)[bus_ix]
+        V_I = get_initial_conditions(sim)[bus_ix + bus_size]
         Vm = sqrt(V_R^2 + V_I^2)
         θ = angle(V_R + V_I * 1im)
         print("Vm ", round(Vm, digits = 4), "\n")
@@ -273,7 +273,12 @@ function show_states_initial_value(sim::Simulation)
         println("====================")
         global_index = global_state_map[name]
         for s in states
-            print(s, " ", round(PSID.get_initial_conditions(sim)[global_index[s]], digits = 4), "\n")
+            print(
+                s,
+                " ",
+                round(get_initial_conditions(sim)[global_index[s]], digits = 4),
+                "\n",
+            )
         end
         println("====================")
     end
@@ -290,7 +295,12 @@ function show_states_initial_value(sim::Simulation)
             global_index = global_state_map[name]
             x0_br = Dict{Symbol, Float64}()
             for (i, s) in enumerate(states)
-                print(s, " ", round(PSID.get_initial_conditions(sim)[global_index[s]], digits = 5), "\n")
+                print(
+                    s,
+                    " ",
+                    round(get_initial_conditions(sim)[global_index[s]], digits = 5),
+                    "\n",
+                )
             end
             println("====================")
         end

From 365c0175aaee87fb6a4147b5545e189edcc74fcb Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:10:57 -0600
Subject: [PATCH 15/22] propagate changes to tests

---
 test/performance/performance_test.jl          |  6 ++----
 test/test_case02_oneDoneQ.jl                  |  6 ++----
 test/test_case03_simple_marconato.jl          | 10 ++++------
 test/test_case04_marconato.jl                 |  6 ++----
 test/test_case05_simple_anderson.jl           |  8 ++------
 test/test_case06_anderson.jl                  |  8 ++------
 test/test_case07_5shaft.jl                    | 11 ++++++-----
 test/test_case08_VirtualSynchMachine.jl       |  2 +-
 test/test_case09_oneDoneQ_inverter.jl         | 11 +++++------
 test/test_case10_staticbranches.jl            |  8 ++++----
 test/test_case11_dynbranches.jl               |  8 ++++----
 test/test_case12_multimachine.jl              |  6 ++++--
 test/test_case13_avrs.jl                      | 11 ++++++-----
 test/test_case14_inverter_ref.jl              |  6 ++++--
 test/test_case15_genrou.jl                    | 10 +++++-----
 test/test_case16_genroe.jl                    | 10 +++++-----
 test/test_case17_genrou_avr.jl                |  5 +++--
 test/test_case18_gensal.jl                    | 10 +++++-----
 test/test_case19_gensae.jl                    | 10 +++++-----
 test/test_case20_ac1a.jl                      | 10 +++++-----
 test/test_case21_gast.jl                      |  6 ++----
 test/test_case22_SteamTurbineGov1.jl          |  6 ++----
 test/test_case23_droopinverter.jl             |  6 ++----
 test/test_case24_gridfollowing.jl             |  4 ++--
 test/test_case25_multimach_dynlines.jl        |  4 ++--
 test/test_case26_SEXS.jl                      | 10 +++++-----
 test/test_case27_source_bus_voltage_change.jl | 16 ++++++----------
 test/test_case28_PeriodicVariableSource.jl    |  6 ++----
 test/test_case29_renewablemodels.jl           | 12 +++++++-----
 test/test_case30_ieeest.jl                    | 10 +++++-----
 test/test_case31_hygov.jl                     |  6 ++----
 test/x_test_sundials.jl                       |  4 ++--
 32 files changed, 115 insertions(+), 137 deletions(-)

diff --git a/test/performance/performance_test.jl b/test/performance/performance_test.jl
index 67154bc38..392054958 100644
--- a/test/performance/performance_test.jl
+++ b/test/performance/performance_test.jl
@@ -21,7 +21,6 @@ try
             m,
             sys,
             pwd(),
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
             console_level = Logging.Error,
         )
@@ -35,11 +34,11 @@ try
         ResidualModel,
         sys,
         pwd(),
-        (0.0, 20.0), #time span
         BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
         console_level = Logging.Error,
     )
-    status = execute!(sim_ida, IDA(), dtmax = 0.01, enable_progress_bar = false)
+    status =
+        execute!(sim_ida, IDA(), (0.0, 20.0), dtmax = 0.01, enable_progress_bar = false)
     if status == PSID.SIMULATION_FINALIZED
         res_ida = read_results(sim_ida)
         solve_time = res_ida.time_log[:timed_solve_time]
@@ -62,7 +61,6 @@ try
         MassMatrixModel,
         sys, #system
         pwd(),
-        (0.0, 20.0), #time span
         BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10"),
         #console_level = Logging.Error,
     ) #Type of Fault
diff --git a/test/test_case02_oneDoneQ.jl b/test/test_case02_oneDoneQ.jl
index b7fa84776..1067d868b 100644
--- a/test/test_case02_oneDoneQ.jl
+++ b/test/test_case02_oneDoneQ.jl
@@ -30,7 +30,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -53,7 +52,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test02_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -89,7 +88,6 @@ end
             MassMatrixModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -112,7 +110,7 @@ end
         @test LinearAlgebra.norm(eigs - test02_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case03_simple_marconato.jl b/test/test_case03_simple_marconato.jl
index 804e8ba25..215d12cfb 100644
--- a/test/test_case03_simple_marconato.jl
+++ b/test/test_case03_simple_marconato.jl
@@ -15,8 +15,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test03.jl"))
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-# time span
-tspan = (0.0, 20.0);
 # Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
     1.0, #change at t = 1.0
@@ -32,7 +30,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -55,7 +53,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test03_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -91,7 +89,7 @@ end
             MassMatrixModel,
             threebus_sys, #system
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -114,7 +112,7 @@ end
         @test LinearAlgebra.norm(eigs - test03_eigvals_psat, Inf) < 5.0
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case04_marconato.jl b/test/test_case04_marconato.jl
index fe441283c..7227f9f4a 100644
--- a/test/test_case04_marconato.jl
+++ b/test/test_case04_marconato.jl
@@ -30,7 +30,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -52,7 +51,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test04_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -88,7 +87,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -110,7 +108,7 @@ end
         @test LinearAlgebra.norm(eigs - test04_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case05_simple_anderson.jl b/test/test_case05_simple_anderson.jl
index 22addbf7f..76565376d 100644
--- a/test/test_case05_simple_anderson.jl
+++ b/test/test_case05_simple_anderson.jl
@@ -15,8 +15,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test05.jl"))
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-# time span
-tspan = (0.0, 200.0);
 # Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
     1.0, #change at t = 1.0
@@ -32,7 +30,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -54,7 +51,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test05_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 200.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -77,7 +74,6 @@ end
             MassMatrixModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -99,7 +95,7 @@ end
         @test LinearAlgebra.norm(eigs - test05_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 200.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case06_anderson.jl b/test/test_case06_anderson.jl
index 9aa6a795c..f92aac84f 100644
--- a/test/test_case06_anderson.jl
+++ b/test/test_case06_anderson.jl
@@ -15,8 +15,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test06.jl"))
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-# time span
-tspan = (0.0, 200.0);
 # Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
     1.0, #change at t = 1.0
@@ -32,7 +30,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -54,7 +51,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test06_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 200.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -77,7 +74,6 @@ end
             MassMatrixModel,
             threebus_sys, #system
             path,
-            (0.0, 20.0), #time span
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -99,7 +95,7 @@ end
         @test LinearAlgebra.norm(eigs - test06_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 200.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case07_5shaft.jl b/test/test_case07_5shaft.jl
index 295adec45..fdf91d23a 100644
--- a/test/test_case07_5shaft.jl
+++ b/test/test_case07_5shaft.jl
@@ -16,7 +16,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test07.jl"))
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-tspan = (0.0, 20.0);
+;
 
 # Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
@@ -33,7 +33,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -55,7 +55,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test07_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.001) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.001) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for angles
@@ -78,7 +78,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -100,7 +100,8 @@ end
         @test LinearAlgebra.norm(eigs - test07_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.001) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.001) ==
+              PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for angles
diff --git a/test/test_case08_VirtualSynchMachine.jl b/test/test_case08_VirtualSynchMachine.jl
index 4105b8034..9a165d958 100644
--- a/test/test_case08_VirtualSynchMachine.jl
+++ b/test/test_case08_VirtualSynchMachine.jl
@@ -54,7 +54,7 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         @test LinearAlgebra.norm(eigs - test08_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case09_oneDoneQ_inverter.jl b/test/test_case09_oneDoneQ_inverter.jl
index ca4c49e48..0c0c50866 100644
--- a/test/test_case09_oneDoneQ_inverter.jl
+++ b/test/test_case09_oneDoneQ_inverter.jl
@@ -18,7 +18,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))
     path = (joinpath(pwd(), "test-09"))
     !isdir(path) && mkdir(path)
     #time span
-    tspan = (0.0, 20.0)
+
     case_inv = collect(PSY.get_components(PSY.DynamicInverter, threebus_sys))[1]
     case_gen = collect(PSY.get_components(PSY.DynamicGenerator, threebus_sys))[1]
 
@@ -31,7 +31,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))
             ResidualModel,
             threebus_sys, # system
             path,
-            tspan,
             [Pref_change, gen_trip],
         )
 
@@ -52,7 +51,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))
         @test LinearAlgebra.norm(eigs - test09_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         #Obtain data for angles
@@ -70,7 +69,7 @@ end
     path = (joinpath(pwd(), "test-09"))
     !isdir(path) && mkdir(path)
     #time span
-    tspan = (0.0, 20.0)
+
     case_inv = collect(PSY.get_components(PSY.DynamicInverter, threebus_sys))[1]
     case_gen = collect(PSY.get_components(PSY.DynamicGenerator, threebus_sys))[1]
 
@@ -83,7 +82,6 @@ end
             MassMatrixModel,
             threebus_sys, # system
             path,
-            tspan,
             [Pref_change, gen_trip],
         )
 
@@ -104,7 +102,8 @@ end
         @test LinearAlgebra.norm(eigs - test09_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02) ==
+              PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         #Obtain data for angles
diff --git a/test/test_case10_staticbranches.jl b/test/test_case10_staticbranches.jl
index 4e12a6fbc..9756b06f2 100644
--- a/test/test_case10_staticbranches.jl
+++ b/test/test_case10_staticbranches.jl
@@ -32,7 +32,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -53,7 +53,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test10_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA()) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), tspan) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for voltages
@@ -79,7 +79,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -100,7 +100,7 @@ end
         @test LinearAlgebra.norm(eigs - test10_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4()) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), tspan) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for voltages
diff --git a/test/test_case11_dynbranches.jl b/test/test_case11_dynbranches.jl
index 356878af2..8f4b6427e 100644
--- a/test/test_case11_dynbranches.jl
+++ b/test/test_case11_dynbranches.jl
@@ -32,7 +32,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
         # Get dictionary of initial conditions
@@ -55,7 +55,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test11_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA()) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), tspan) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for voltages
@@ -75,7 +75,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -96,7 +96,7 @@ end
         @test LinearAlgebra.norm(eigs - test11_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas5()) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas5(), tspan) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for voltages
diff --git a/test/test_case12_multimachine.jl b/test/test_case12_multimachine.jl
index 25cb32016..1009c67dd 100644
--- a/test/test_case12_multimachine.jl
+++ b/test/test_case12_multimachine.jl
@@ -33,7 +33,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -58,6 +58,7 @@ Ybus_change = NetworkSwitch(
         @test execute!(
             sim, #simulation structure
             IDA(),#Sundials DAE Solver
+            tspan,
             dtmax = 0.02, #keywords arguments
         ) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
@@ -77,7 +78,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -102,6 +103,7 @@ end
         @test execute!(
             sim, #simulation structure
             Rodas4(),#Sundials DAE Solver
+            tspan,
             dtmax = 0.02, #keywords arguments
         ) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
diff --git a/test/test_case13_avrs.jl b/test/test_case13_avrs.jl
index c15296bb2..e02bd58ef 100644
--- a/test/test_case13_avrs.jl
+++ b/test/test_case13_avrs.jl
@@ -15,7 +15,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test13.jl"))
 ##################################################
 
 # Time span
-tspan = (0.0, 20.0)
+
 
 #Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
@@ -32,7 +32,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -54,7 +54,7 @@ Ybus_change = NetworkSwitch(
         @test LinearAlgebra.norm(eigs - test13_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         #Obtain data for angles
@@ -74,7 +74,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -96,7 +96,8 @@ end
         @test LinearAlgebra.norm(eigs - test13_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.02) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02) ==
+              PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         #Obtain data for angles
diff --git a/test/test_case14_inverter_ref.jl b/test/test_case14_inverter_ref.jl
index ec9f329f3..257000fec 100644
--- a/test/test_case14_inverter_ref.jl
+++ b/test/test_case14_inverter_ref.jl
@@ -33,7 +33,7 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -58,6 +58,7 @@ Ybus_change = NetworkSwitch(
         @test execute!(
             sim, #simulation structure
             IDA(),#Sundials DAE Solver
+            tspan,
             dtmax = 0.001, #keywords arguments
         ) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
@@ -78,7 +79,7 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -103,6 +104,7 @@ end
         @test execute!(
             sim, #simulation structure
             Rodas4(),#Sundials DAE Solver
+            tspan,
             dtmax = 0.001, #keywords arguments
         ) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
diff --git a/test/test_case15_genrou.jl b/test/test_case15_genrou.jl
index 52522b8a0..2dd4d4a7b 100644
--- a/test/test_case15_genrou.jl
+++ b/test/test_case15_genrou.jl
@@ -30,7 +30,7 @@ init_conditions =
 eigs_values = [test15_eigvals, test15_eigvals_no_sat, test15_eigvals_high_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROU/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-genrou"))
@@ -43,7 +43,7 @@ function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -65,7 +65,7 @@ function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -104,7 +104,7 @@ function test_genrou_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -126,7 +126,7 @@ function test_genrou_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case16_genroe.jl b/test/test_case16_genroe.jl
index 582d67e8c..5e637f7d6 100644
--- a/test/test_case16_genroe.jl
+++ b/test/test_case16_genroe.jl
@@ -27,7 +27,7 @@ init_conditions = [test_psse_genroe_init, test_psse_genroe_high_sat_init]
 eigs_values = [test16_eigvals, test16_eigvals_high_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-genrou"))
@@ -40,7 +40,7 @@ function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -62,7 +62,7 @@ function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -101,7 +101,7 @@ function test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -123,7 +123,7 @@ function test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case17_genrou_avr.jl b/test/test_case17_genrou_avr.jl
index c945e5e23..5e942b840 100644
--- a/test/test_case17_genrou_avr.jl
+++ b/test/test_case17_genrou_avr.jl
@@ -45,7 +45,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test17.jl"))
         @test LinearAlgebra.norm(eigs - test17_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.01) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.01) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for angles
@@ -87,7 +87,8 @@ end
         @test LinearAlgebra.norm(eigs - test17_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.01) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.01) ==
+              PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for angles
diff --git a/test/test_case18_gensal.jl b/test/test_case18_gensal.jl
index a9c9ca94d..9df460ac3 100644
--- a/test/test_case18_gensal.jl
+++ b/test/test_case18_gensal.jl
@@ -37,7 +37,7 @@ init_conditions = [
 eigs_values = [test18_eigvals]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENSAL/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-gensal"))
@@ -50,7 +50,7 @@ function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -72,7 +72,7 @@ function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -111,7 +111,7 @@ function test_gensal_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -133,7 +133,7 @@ function test_gensal_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case19_gensae.jl b/test/test_case19_gensae.jl
index 933cc9f99..7c0b4f073 100644
--- a/test/test_case19_gensae.jl
+++ b/test/test_case19_gensae.jl
@@ -33,7 +33,7 @@ init_conditions = [
 eigs_values = [test19_eigvals]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENSAE/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-gensae"))
@@ -46,7 +46,7 @@ function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -68,7 +68,7 @@ function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -107,7 +107,7 @@ function test_gensae_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -129,7 +129,7 @@ function test_gensae_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case20_ac1a.jl b/test/test_case20_ac1a.jl
index fe041de33..8d9c8d1b4 100644
--- a/test/test_case20_ac1a.jl
+++ b/test/test_case20_ac1a.jl
@@ -28,7 +28,7 @@ init_conditions = [test_psse_ac1a_init, test_psse_ac1a_sat_init]
 eigs_values = [test20_eigvals, test20_eigvals_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/AC1A/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-ac1a"))
@@ -41,7 +41,7 @@ function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -63,7 +63,7 @@ function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -101,7 +101,7 @@ function test_ac1a_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -123,7 +123,7 @@ function test_ac1a_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case21_gast.jl b/test/test_case21_gast.jl
index bd6bcc7bc..b8a3fd86d 100644
--- a/test/test_case21_gast.jl
+++ b/test/test_case21_gast.jl
@@ -23,7 +23,6 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/GAST/GAST_TEST.csv")
             ResidualModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -45,7 +44,7 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/GAST/GAST_TEST.csv")
         @test LinearAlgebra.norm(eigs - test21_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -78,7 +77,6 @@ end
             MassMatrixModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -100,7 +98,7 @@ end
         @test LinearAlgebra.norm(eigs - test21_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case22_SteamTurbineGov1.jl b/test/test_case22_SteamTurbineGov1.jl
index 0695118cb..bb5e0003a 100644
--- a/test/test_case22_SteamTurbineGov1.jl
+++ b/test/test_case22_SteamTurbineGov1.jl
@@ -23,7 +23,6 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/TGOV1/TEST_TGOV1.csv")
             ResidualModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -45,7 +44,7 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/TGOV1/TEST_TGOV1.csv")
         @test LinearAlgebra.norm(eigs - test22_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -78,7 +77,6 @@ end
             MassMatrixModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -100,7 +98,7 @@ end
         @test LinearAlgebra.norm(eigs - test22_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case23_droopinverter.jl b/test/test_case23_droopinverter.jl
index 7734aac5e..dc16242c5 100644
--- a/test/test_case23_droopinverter.jl
+++ b/test/test_case23_droopinverter.jl
@@ -36,7 +36,6 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
             ResidualModel,
             omib_sys, # system
             path,
-            tspan,
             Pref_change,
         )
 
@@ -57,7 +56,7 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         @test LinearAlgebra.norm(eigs - test23_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), tspan, dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -90,7 +89,6 @@ end
             MassMatrixModel,
             omib_sys, # system
             path,
-            tspan,
             Pref_change,
         )
 
@@ -111,7 +109,7 @@ end
         @test LinearAlgebra.norm(eigs - test23_eigvals) < 1e-3
 
         #Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), tspan, dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case24_gridfollowing.jl b/test/test_case24_gridfollowing.jl
index 87df0d5d1..89f1f0c2b 100644
--- a/test/test_case24_gridfollowing.jl
+++ b/test/test_case24_gridfollowing.jl
@@ -57,7 +57,7 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3
 
         #Solve problem in equilibrium
-        @test execute!(sim, Sundials.IDA(), dtmax = 0.001, saveat = 0.005) ==
+        @test execute!(sim, Sundials.IDA(), (0.0, 2.0), dtmax = 0.001, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -118,7 +118,7 @@ end
         @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3
 
         #Solve problem in equilibrium
-        @test execute!(sim, Rodas4(), dtmax = 0.001, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 2.0), dtmax = 0.001, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case25_multimach_dynlines.jl b/test/test_case25_multimach_dynlines.jl
index cd7103042..e6109de2d 100644
--- a/test/test_case25_multimach_dynlines.jl
+++ b/test/test_case25_multimach_dynlines.jl
@@ -48,7 +48,7 @@ Pref_change = ControlReferenceChange(1.0, gen2, :P_ref, 0.9);
         @test small_sig.stable
 
         # Solve problem in equilibrium
-        @test execute!(sim, Sundials.IDA(), dtmax = 0.01, saveat = 0.01) ==
+        @test execute!(sim, Sundials.IDA(), (0.0, 40.0), dtmax = 0.01, saveat = 0.01) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -91,7 +91,7 @@ end
         @test small_sig.stable
 
         # Solve problem in equilibrium
-        @test execute!(sim, Rodas4(), dtmax = 0.01, saveat = 0.01) ==
+        @test execute!(sim, Rodas4(), (0.0, 40.0), dtmax = 0.01, saveat = 0.01) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case26_SEXS.jl b/test/test_case26_SEXS.jl
index c2ba43157..7275d450f 100644
--- a/test/test_case26_SEXS.jl
+++ b/test/test_case26_SEXS.jl
@@ -28,7 +28,7 @@ init_conditions = test26_x0_init
 eigs_values = [test26_eigvals, test26_eigvals_noTE]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/SEXS/ThreeBusMulti.raw")
-tspan = (0.0, 20.0)
+
 
 function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-sexs"))
@@ -41,7 +41,7 @@ function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -63,7 +63,7 @@ function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -103,7 +103,7 @@ function test_sexs_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -125,7 +125,7 @@ function test_sexs_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case27_source_bus_voltage_change.jl b/test/test_case27_source_bus_voltage_change.jl
index 8ef8ac647..57fd2d5ba 100644
--- a/test/test_case27_source_bus_voltage_change.jl
+++ b/test/test_case27_source_bus_voltage_change.jl
@@ -18,7 +18,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))
 ####### Changing magnitude of votlage at source bus #########
 
 # time span
-tspan = (0.0, 20.0);
+;
 case_source = collect(PSY.get_components(PSY.Source, threebus_sys))[1]
 
 # Define Fault using Callbacks
@@ -33,7 +33,6 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.V_source_index,
             ResidualModel,
             threebus_sys, # system
             path,
-            tspan,
             V_source_change,
         )
 
@@ -46,7 +45,7 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.V_source_index,
         @test (diff[1] < 1e-3)
 
         # Solve problem
-        execute!(sim, IDA(), dtmax = 0.02)
+        execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02)
         results = read_results(sim)
 
         # Obtain data for angles
@@ -66,7 +65,6 @@ end
             MassMatrixModel,
             threebus_sys, # system
             path,
-            tspan,
             V_source_change,
         )
 
@@ -79,7 +77,7 @@ end
         @test (diff[1] < 1e-3)
 
         # Solve problem
-        execute!(sim, Rodas4(), dtmax = 0.02)
+        execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02)
         results = read_results(sim)
 
         # Obtain data for angles
@@ -93,7 +91,7 @@ end
 ####### Changing angle of voltage at source bus #########
 
 #time span
-tspan = (0.0, 20.0);
+;
 case_source = collect(PSY.get_components(PSY.Source, threebus_sys))[1]
 
 #Define Fault using Callbacks
@@ -108,7 +106,6 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.θ_source_index,
             ResidualModel,
             threebus_sys, # system
             path,
-            tspan,
             V_source_change,
         )
 
@@ -121,7 +118,7 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.θ_source_index,
         @test (diff[1] < 1e-3)
 
         # Solve problem
-        execute!(sim, IDA(), dtmax = 0.02)
+        execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02)
         results = read_results(sim)
 
         # Obtain data for angles
@@ -140,7 +137,6 @@ end
             MassMatrixModel,
             threebus_sys, # system
             path,
-            tspan,
             V_source_change,
         )
 
@@ -153,7 +149,7 @@ end
         @test (diff[1] < 1e-3)
 
         # Solve problem
-        execute!(sim, Rodas4(), dtmax = 0.02)
+        execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02)
         results = read_results(sim)
 
         # Obtain data for angles
diff --git a/test/test_case28_PeriodicVariableSource.jl b/test/test_case28_PeriodicVariableSource.jl
index d2d130bbb..9bf327a36 100644
--- a/test/test_case28_PeriodicVariableSource.jl
+++ b/test/test_case28_PeriodicVariableSource.jl
@@ -27,11 +27,10 @@ tsteps = tspan[1]:step:tspan[2]
             ResidualModel,
             sys, # system
             path,
-            tspan,
         )
 
         #Solve problem
-        @test execute!(sim, IDA(), saveat = tsteps) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), tspan, saveat = tsteps) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
         pvs = collect(get_components(PeriodicVariableSource, sys))[1]
 
@@ -68,11 +67,10 @@ end
             MassMatrixModel,
             sys, # system
             path,
-            tspan,
         )
 
         #Solve problem
-        @test execute!(sim, Rodas4(), saveat = tsteps) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, Rodas4(), tspan, saveat = tsteps) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
         pvs = collect(get_components(PeriodicVariableSource, sys))[1]
 
diff --git a/test/test_case29_renewablemodels.jl b/test/test_case29_renewablemodels.jl
index 3f0f84c9d..4097fac0e 100644
--- a/test/test_case29_renewablemodels.jl
+++ b/test/test_case29_renewablemodels.jl
@@ -67,6 +67,7 @@ function test_renA_implicit(csv_file, init_cond, eigs_value, F_Flag)
         @test execute!(
             sim,
             IDA(),
+            tspan,
             dtmax = 0.005,
             saveat = 0.005,
             abstol = 1e-9,
@@ -137,6 +138,7 @@ function test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_Flag)
         @test execute!(
             sim,
             Rodas4(),
+            tspan,
             dtmax = 0.005,
             saveat = 0.005,
             abstol = 1e-6,
@@ -200,7 +202,7 @@ end
 Case 29: Test with dyr files
 This case study a three bus system with 1 Generic Renewable Model (REPCA, REECB, REGCA) and an infinite source.
 The fault drop the connection between buses 1 and 3, eliminating the direct connection between the load at bus 1
-and the generator located in bus 3. The infinite generator is located at bus 2. 
+and the generator located in bus 3. The infinite generator is located at bus 2.
 """
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/ThreeBusRenewable.raw")
@@ -251,7 +253,7 @@ function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
         ) #Type of Fault
 
@@ -273,7 +275,7 @@ function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -314,7 +316,7 @@ function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tsp
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
         ) #Type of Fault
 
@@ -336,7 +338,7 @@ function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tsp
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case30_ieeest.jl b/test/test_case30_ieeest.jl
index 8970ce73f..0554d5c13 100644
--- a/test/test_case30_ieeest.jl
+++ b/test/test_case30_ieeest.jl
@@ -25,7 +25,7 @@ init_conditions = [test_psse_ieeest_no_filt_init, test_psse_ieeest_with_filt_ini
 
 eigs_values = [test30_eigvals_no_filt, test30_eigvals_with_filt]
 
-tspan = (0.0, 20.0)
+
 
 function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-ieeest"))
@@ -38,7 +38,7 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -60,7 +60,7 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -105,7 +105,7 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-            tspan, #time span
+
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -127,7 +127,7 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case31_hygov.jl b/test/test_case31_hygov.jl
index 631dac679..8b616a210 100644
--- a/test/test_case31_hygov.jl
+++ b/test/test_case31_hygov.jl
@@ -23,7 +23,6 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/HYGOV_RESULTS.csv")
             ResidualModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -45,7 +44,7 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/HYGOV_RESULTS.csv")
         @test LinearAlgebra.norm(eigs - test31_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -84,7 +83,6 @@ end
             MassMatrixModel,
             sys, #system
             path,
-            (0.0, 20.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         )
 
@@ -106,7 +104,7 @@ end
         @test LinearAlgebra.norm(eigs - test31_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/x_test_sundials.jl b/test/x_test_sundials.jl
index 5b46200ef..b29019caf 100644
--- a/test/x_test_sundials.jl
+++ b/test/x_test_sundials.jl
@@ -44,7 +44,7 @@ function test_sundials(solver)
             ResidualModel,
             path,
             threebus_sys, #system
-            tspan, #time span
+
             Ybus_change, #Type of Fault
         )
 
@@ -52,7 +52,7 @@ function test_sundials(solver)
         small_sig = small_signal_analysis(sim)
 
         #Solve problem
-        @info "$(solver)" @time execute!(sim, IDA(linear_solver = solver);)
+        @info "$(solver)" @time execute!(sim, IDA(linear_solver = solver), tspan)
 
         #Obtain data for voltages
         series = get_voltage_magnitude_series(results, 102)

From 46bba00ebb2037e09b19c56177f6fc1885d67985 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:19:21 -0600
Subject: [PATCH 16/22] fix more tests

---
 test/test_case03_simple_marconato.jl    | 2 --
 test/test_case07_5shaft.jl              | 2 --
 test/test_case08_VirtualSynchMachine.jl | 6 ++----
 test/test_case10_staticbranches.jl      | 2 --
 test/test_case11_dynbranches.jl         | 2 --
 test/test_case12_multimachine.jl        | 2 --
 test/test_case13_avrs.jl                | 3 ---
 test/test_case14_inverter_ref.jl        | 2 --
 test/test_case15_genrou.jl              | 3 ---
 test/test_case16_genroe.jl              | 3 ---
 test/test_case17_genrou_avr.jl          | 6 ++----
 test/test_case18_gensal.jl              | 3 ---
 test/test_case19_gensae.jl              | 3 ---
 test/test_case20_ac1a.jl                | 3 ---
 test/test_case25_multimach_dynlines.jl  | 4 ++--
 test/test_case26_SEXS.jl                | 3 ---
 test/test_case29_renewablemodels.jl     | 6 ++----
 test/test_case30_ieeest.jl              | 4 ----
 test/x_test_sundials.jl                 | 1 -
 19 files changed, 8 insertions(+), 52 deletions(-)

diff --git a/test/test_case03_simple_marconato.jl b/test/test_case03_simple_marconato.jl
index 215d12cfb..6ef944923 100644
--- a/test/test_case03_simple_marconato.jl
+++ b/test/test_case03_simple_marconato.jl
@@ -30,7 +30,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system
             path,
-
             Ybus_change, #Type of Fault
         )
 
@@ -89,7 +88,6 @@ end
             MassMatrixModel,
             threebus_sys, #system
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case07_5shaft.jl b/test/test_case07_5shaft.jl
index fdf91d23a..f752760cf 100644
--- a/test/test_case07_5shaft.jl
+++ b/test/test_case07_5shaft.jl
@@ -33,7 +33,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         ) #initial guess
 
@@ -78,7 +77,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         ) #initial guess
 
diff --git a/test/test_case08_VirtualSynchMachine.jl b/test/test_case08_VirtualSynchMachine.jl
index 9a165d958..110cc4b83 100644
--- a/test/test_case08_VirtualSynchMachine.jl
+++ b/test/test_case08_VirtualSynchMachine.jl
@@ -33,7 +33,6 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
             ResidualModel,
             omib_sys, # system
             path,
-            (0.0, 4.0),
             Pref_change,
         )
 
@@ -54,7 +53,7 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         @test LinearAlgebra.norm(eigs - test08_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), (0.0, 4.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -90,7 +89,6 @@ end
             MassMatrixModel,
             omib_sys, # system
             path,
-            (0.0, 4.0),
             Pref_change,
         )
 
@@ -111,7 +109,7 @@ end
         @test LinearAlgebra.norm(eigs - test08_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas5(), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas5(), (0.0, 4.0), dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case10_staticbranches.jl b/test/test_case10_staticbranches.jl
index 9756b06f2..b5753c9e1 100644
--- a/test/test_case10_staticbranches.jl
+++ b/test/test_case10_staticbranches.jl
@@ -32,7 +32,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
@@ -79,7 +78,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case11_dynbranches.jl b/test/test_case11_dynbranches.jl
index 8f4b6427e..659f7d09a 100644
--- a/test/test_case11_dynbranches.jl
+++ b/test/test_case11_dynbranches.jl
@@ -32,7 +32,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
         # Get dictionary of initial conditions
@@ -75,7 +74,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case12_multimachine.jl b/test/test_case12_multimachine.jl
index 1009c67dd..386c4ba02 100644
--- a/test/test_case12_multimachine.jl
+++ b/test/test_case12_multimachine.jl
@@ -33,7 +33,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
@@ -78,7 +77,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case13_avrs.jl b/test/test_case13_avrs.jl
index e02bd58ef..a79c9cc14 100644
--- a/test/test_case13_avrs.jl
+++ b/test/test_case13_avrs.jl
@@ -16,7 +16,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test13.jl"))
 
 # Time span
 
-
 #Define Fault: Change of YBus
 Ybus_change = NetworkSwitch(
     1.0, #change at t = 1.0
@@ -32,7 +31,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
@@ -74,7 +72,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case14_inverter_ref.jl b/test/test_case14_inverter_ref.jl
index 257000fec..20bab6524 100644
--- a/test/test_case14_inverter_ref.jl
+++ b/test/test_case14_inverter_ref.jl
@@ -33,7 +33,6 @@ Ybus_change = NetworkSwitch(
             ResidualModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
@@ -79,7 +78,6 @@ end
             MassMatrixModel,
             threebus_sys, #system,
             path,
-
             Ybus_change, #Type of Fault
         )
 
diff --git a/test/test_case15_genrou.jl b/test/test_case15_genrou.jl
index 2dd4d4a7b..d7f35f95b 100644
--- a/test/test_case15_genrou.jl
+++ b/test/test_case15_genrou.jl
@@ -31,7 +31,6 @@ eigs_values = [test15_eigvals, test15_eigvals_no_sat, test15_eigvals_high_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROU/ThreeBusMulti.raw")
 
-
 function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-genrou"))
     !isdir(path) && mkdir(path)
@@ -43,7 +42,6 @@ function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -104,7 +102,6 @@ function test_genrou_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case16_genroe.jl b/test/test_case16_genroe.jl
index 5e637f7d6..a72550e22 100644
--- a/test/test_case16_genroe.jl
+++ b/test/test_case16_genroe.jl
@@ -28,7 +28,6 @@ eigs_values = [test16_eigvals, test16_eigvals_high_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENROE/ThreeBusMulti.raw")
 
-
 function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-genrou"))
     !isdir(path) && mkdir(path)
@@ -40,7 +39,6 @@ function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -101,7 +99,6 @@ function test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case17_genrou_avr.jl b/test/test_case17_genrou_avr.jl
index 5e942b840..7536bad39 100644
--- a/test/test_case17_genrou_avr.jl
+++ b/test/test_case17_genrou_avr.jl
@@ -23,7 +23,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test17.jl"))
             ResidualModel,
             sys, #system
             path,
-            (0.0, 30.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault,
         ) #Type of Fault
 
@@ -45,7 +44,7 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test17.jl"))
         @test LinearAlgebra.norm(eigs - test17_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.01) == PSID.SIMULATION_FINALIZED
+        @test execute!(sim, IDA(), (0.0, 30.0), dtmax = 0.01) == PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
         # Obtain data for angles
@@ -65,7 +64,6 @@ end
             MassMatrixModel,
             sys, #system
             path,
-            (0.0, 30.0), #time span
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault,
         ) #Type of Fault
 
@@ -87,7 +85,7 @@ end
         @test LinearAlgebra.norm(eigs - test17_eigvals) < 1e-3
 
         # Solve problem
-        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.01) ==
+        @test execute!(sim, Rodas4(), (0.0, 30.0), dtmax = 0.01) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
diff --git a/test/test_case18_gensal.jl b/test/test_case18_gensal.jl
index 9df460ac3..2e5507f0b 100644
--- a/test/test_case18_gensal.jl
+++ b/test/test_case18_gensal.jl
@@ -38,7 +38,6 @@ eigs_values = [test18_eigvals]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENSAL/ThreeBusMulti.raw")
 
-
 function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-gensal"))
     !isdir(path) && mkdir(path)
@@ -50,7 +49,6 @@ function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -111,7 +109,6 @@ function test_gensal_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case19_gensae.jl b/test/test_case19_gensae.jl
index 7c0b4f073..46a6af05c 100644
--- a/test/test_case19_gensae.jl
+++ b/test/test_case19_gensae.jl
@@ -34,7 +34,6 @@ eigs_values = [test19_eigvals]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/GENSAE/ThreeBusMulti.raw")
 
-
 function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-gensae"))
     !isdir(path) && mkdir(path)
@@ -46,7 +45,6 @@ function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -107,7 +105,6 @@ function test_gensae_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case20_ac1a.jl b/test/test_case20_ac1a.jl
index 8d9c8d1b4..635fdbd24 100644
--- a/test/test_case20_ac1a.jl
+++ b/test/test_case20_ac1a.jl
@@ -29,7 +29,6 @@ eigs_values = [test20_eigvals, test20_eigvals_sat]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/AC1A/ThreeBusMulti.raw")
 
-
 function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-ac1a"))
     !isdir(path) && mkdir(path)
@@ -41,7 +40,6 @@ function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -101,7 +99,6 @@ function test_ac1a_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case25_multimach_dynlines.jl b/test/test_case25_multimach_dynlines.jl
index e6109de2d..11682ce7d 100644
--- a/test/test_case25_multimach_dynlines.jl
+++ b/test/test_case25_multimach_dynlines.jl
@@ -33,7 +33,7 @@ Pref_change = ControlReferenceChange(1.0, gen2, :P_ref, 0.9);
     !isdir(path) && mkdir(path)
     try
         # Define Simulation Problem
-        sim = Simulation!(ResidualModel, sys, path, tspan, Pref_change)
+        sim = Simulation!(ResidualModel, sys, path, Pref_change)
 
         # Test Initial Condition
         diff = [0.0]
@@ -76,7 +76,7 @@ end
     !isdir(path) && mkdir(path)
     try
         # Define Simulation Problem
-        sim = Simulation!(MassMatrixModel, sys, path, tspan, Pref_change)
+        sim = Simulation!(MassMatrixModel, sys, path, Pref_change)
 
         # Test Initial Condition
         diff = [0.0]
diff --git a/test/test_case26_SEXS.jl b/test/test_case26_SEXS.jl
index 7275d450f..8a573c57d 100644
--- a/test/test_case26_SEXS.jl
+++ b/test/test_case26_SEXS.jl
@@ -29,7 +29,6 @@ eigs_values = [test26_eigvals, test26_eigvals_noTE]
 
 raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/SEXS/ThreeBusMulti.raw")
 
-
 function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-sexs"))
     !isdir(path) && mkdir(path)
@@ -41,7 +40,6 @@ function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -103,7 +101,6 @@ function test_sexs_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case29_renewablemodels.jl b/test/test_case29_renewablemodels.jl
index 4097fac0e..65c6a4813 100644
--- a/test/test_case29_renewablemodels.jl
+++ b/test/test_case29_renewablemodels.jl
@@ -44,7 +44,7 @@ function test_renA_implicit(csv_file, init_cond, eigs_value, F_Flag)
 
         Ybus_change = BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2")
 
-        sim = Simulation(ResidualModel, sys, path, tspan, Ybus_change)
+        sim = Simulation(ResidualModel, sys, path, Ybus_change)
 
         # Test Initial Condition
         diff = [0.0]
@@ -115,7 +115,7 @@ function test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_Flag)
 
         Ybus_change = BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2")
 
-        sim = Simulation(MassMatrixModel, sys, path, tspan, Ybus_change)
+        sim = Simulation(MassMatrixModel, sys, path, Ybus_change)
 
         # Test Initial Condition
         diff = [0.0]
@@ -253,7 +253,6 @@ function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
         ) #Type of Fault
 
@@ -316,7 +315,6 @@ function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tsp
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/test_case30_ieeest.jl b/test/test_case30_ieeest.jl
index 0554d5c13..193cee311 100644
--- a/test/test_case30_ieeest.jl
+++ b/test/test_case30_ieeest.jl
@@ -25,8 +25,6 @@ init_conditions = [test_psse_ieeest_no_filt_init, test_psse_ieeest_with_filt_ini
 
 eigs_values = [test30_eigvals_no_filt, test30_eigvals_with_filt]
 
-
-
 function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
     path = (joinpath(pwd(), "test-psse-ieeest"))
     !isdir(path) && mkdir(path)
@@ -38,7 +36,6 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
             ResidualModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
@@ -105,7 +102,6 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
             MassMatrixModel,
             sys, #system
             path,
-
             BranchTrip(1.0, Line, "BUS 1-BUS 2-i_1"), #Type of Fault
         ) #Type of Fault
 
diff --git a/test/x_test_sundials.jl b/test/x_test_sundials.jl
index b29019caf..998cd64ae 100644
--- a/test/x_test_sundials.jl
+++ b/test/x_test_sundials.jl
@@ -44,7 +44,6 @@ function test_sundials(solver)
             ResidualModel,
             path,
             threebus_sys, #system
-
             Ybus_change, #Type of Fault
         )
 

From 270cc676df610efa813c052bacc2e12b24bf6bd0 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 15:29:30 -0600
Subject: [PATCH 17/22] fix test

---
 test/test_case24_gridfollowing.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/test/test_case24_gridfollowing.jl b/test/test_case24_gridfollowing.jl
index 89f1f0c2b..02bd1b965 100644
--- a/test/test_case24_gridfollowing.jl
+++ b/test/test_case24_gridfollowing.jl
@@ -17,9 +17,6 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test24.jl"))
 ############### SOLVE PROBLEM ####################
 ##################################################
 
-# time span
-tspan = (0.0, 2.0);
-
 # PSCAD benchmark data
 csv_file = joinpath(TEST_FILES_DIR, "benchmarks/pscad/Test24/Test24_p.csv")
 t_offset = 9.0
@@ -36,7 +33,6 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
             ResidualModel,
             omib_sys, # system
             path,
-            tspan,
             Pref_change,
         )
 

From b8c17fe4fe5ee56a9c2420c4710d4a15c282680e Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 16:01:35 -0600
Subject: [PATCH 18/22] more tests fixes

---
 test/test_case24_gridfollowing.jl   | 1 -
 test/test_case29_renewablemodels.jl | 4 ++--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/test/test_case24_gridfollowing.jl b/test/test_case24_gridfollowing.jl
index 02bd1b965..0bf47ae09 100644
--- a/test/test_case24_gridfollowing.jl
+++ b/test/test_case24_gridfollowing.jl
@@ -93,7 +93,6 @@ end
             MassMatrixModel,
             omib_sys, # system
             path,
-            tspan,
             Pref_change,
         )
 
diff --git a/test/test_case29_renewablemodels.jl b/test/test_case29_renewablemodels.jl
index 65c6a4813..9ca21bf95 100644
--- a/test/test_case29_renewablemodels.jl
+++ b/test/test_case29_renewablemodels.jl
@@ -274,7 +274,7 @@ function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2
 
         # Solve problem
-        @test execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, IDA(), tspan, dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 
@@ -336,7 +336,7 @@ function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tsp
         @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2
 
         # Solve problem
-        @test execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.005, saveat = 0.005) ==
+        @test execute!(sim, Rodas4(), tspan, dtmax = 0.005, saveat = 0.005) ==
               PSID.SIMULATION_FINALIZED
         results = read_results(sim)
 

From f853c8fca58ba41f3502026cfd064b1aa8b24d0b Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 16:22:42 -0600
Subject: [PATCH 19/22] reorder initialization code

---
 src/base/simulation.jl                | 18 +++++-------------
 src/base/simulation_initialization.jl | 24 ++++++++++++++++--------
 2 files changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/base/simulation.jl b/src/base/simulation.jl
index 0d68967f3..9a5b1934d 100644
--- a/src/base/simulation.jl
+++ b/src/base/simulation.jl
@@ -207,22 +207,14 @@ function _build_inputs!(
     return
 end
 
-function _get_flat_start(inputs::SimulationInputs)
-    bus_count = get_bus_count(inputs)
-    var_count = get_variable_count(inputs)
-    initial_conditions = zeros(var_count)
-    initial_conditions[1:bus_count] .= 1.0
-    return initial_conditions
-end
-
 function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel}
     simulation_inputs = get_simulation_inputs(sim)
     x0_init = deepcopy(sim.user_provided_initial_conditions)
     if isempty(x0_init) && sim.initialized
         @warn "Initial Conditions set to flat start"
-        x0_init = _get_flat_start(simulation_inputs)
+        x0_init = get_flat_start(simulation_inputs)
     elseif isempty(x0_init) && !sim.initialized
-        x0_init = _get_flat_start(simulation_inputs)
+        x0_init = get_flat_start(simulation_inputs)
     elseif !isempty(x0_init) && sim.initialized
         if length(x0_init) != get_variable_count(simulation_inputs)
             throw(
@@ -234,7 +226,7 @@ function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel
     elseif !isempty(x0_init) && !sim.initialized
         @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
         if length(x0_init) != get_variable_count(simulation_inputs)
-            x0_init = _get_flat_start(simulation_inputs)
+            x0_init = get_flat_start(simulation_inputs)
         end
     else
         @assert false
@@ -295,10 +287,10 @@ function _get_model_function(
 )
     return SciMLBase.DAEFunction{true}(
         model;
-        # Currently commented for Sundials compatibility
+        # Currently commented for Sundials.jl compatibility
         #jac = jacobian,
         tgrad = (dT, u, p, t) -> dT .= false,
-        jac_prototype = jacobian.Jv,
+        #jac_prototype = jacobian.Jv,
     )
 end
 
diff --git a/src/base/simulation_initialization.jl b/src/base/simulation_initialization.jl
index 4a7216bb5..68074d8a6 100644
--- a/src/base/simulation_initialization.jl
+++ b/src/base/simulation_initialization.jl
@@ -1,4 +1,12 @@
-function _power_flow_solution!(
+function get_flat_start(inputs::SimulationInputs)
+    bus_count = get_bus_count(inputs)
+    var_count = get_variable_count(inputs)
+    initial_conditions = zeros(var_count)
+    initial_conditions[1:bus_count] .= 1.0
+    return initial_conditions
+end
+
+function power_flow_solution!(
     initial_guess::Vector{Float64},
     sys::PSY.System,
     inputs::SimulationInputs,
@@ -20,7 +28,7 @@ function _power_flow_solution!(
     return BUILD_INCOMPLETE
 end
 
-function _initialize_static_injection!(inputs::SimulationInputs)
+function initialize_static_injection!(inputs::SimulationInputs)
     @debug "Updating Source internal voltage magnitude and angle"
     static_injection_devices = get_static_injectors(inputs)
     if !isempty(static_injection_devices)
@@ -37,7 +45,7 @@ function _initialize_static_injection!(inputs::SimulationInputs)
     return BUILD_INCOMPLETE
 end
 
-function _initialize_dynamic_injection!(
+function initialize_dynamic_injection!(
     initial_guess::Vector{Float64},
     inputs::SimulationInputs,
     system::PSY.System,
@@ -67,7 +75,7 @@ function _initialize_dynamic_injection!(
     return BUILD_INCOMPLETE
 end
 
-function _initialize_dynamic_branches!(
+function initialize_dynamic_branches!(
     initial_guess::Vector{Float64},
     inputs::SimulationInputs,
 )
@@ -130,17 +138,17 @@ function _calculate_initial_guess!(x0_init::Vector{Float64}, sim::Simulation)
     while sim.status == BUILD_INCOMPLETE
         @debug "Start state intialization routine"
         TimerOutputs.@timeit BUILD_TIMER "Power Flow solution" begin
-            sim.status = _power_flow_solution!(x0_init, get_system(sim), inputs)
+            sim.status = power_flow_solution!(x0_init, get_system(sim), inputs)
         end
         TimerOutputs.@timeit BUILD_TIMER "Initialize Static Injectors" begin
-            sim.status = _initialize_static_injection!(inputs)
+            sim.status = initialize_static_injection!(inputs)
         end
         TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Injectors" begin
-            sim.status = _initialize_dynamic_injection!(x0_init, inputs, get_system(sim))
+            sim.status = initialize_dynamic_injection!(x0_init, inputs, get_system(sim))
         end
         if has_dyn_lines(inputs)
             TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Branches" begin
-                sim.status = _initialize_dynamic_branches!(x0_init, inputs)
+                sim.status = initialize_dynamic_branches!(x0_init, inputs)
             end
         else
             @debug "No Dynamic Branches in the system"

From 74f3e16a7afb7c38194847c8f3d0909b402ad02d Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 16:46:04 -0600
Subject: [PATCH 20/22] add functions to get jacobian

---
 src/PowerSimulationsDynamics.jl | 17 +++++++------
 src/base/jacobian.jl            | 43 +++++++++++++++++++++++++--------
 test/test_base.jl               |  6 +++++
 3 files changed, 48 insertions(+), 18 deletions(-)

diff --git a/src/PowerSimulationsDynamics.jl b/src/PowerSimulationsDynamics.jl
index 25cf7aa2e..d881a67d2 100644
--- a/src/PowerSimulationsDynamics.jl
+++ b/src/PowerSimulationsDynamics.jl
@@ -30,17 +30,18 @@ export LoadChange
 # export BusTrip
 
 # Export for routines
-export small_signal_analysis
-export get_state_series
-export get_voltage_magnitude_series
-export get_voltage_angle_series
-export show_states_initial_value
-export read_initial_conditions
-export get_real_current_series
-export get_imaginary_current_series
 export get_activepower_series
+export get_jacobian
+export get_imaginary_current_series
 export get_reactivepower_series
+export get_real_current_series
 export get_setpoints
+export get_state_series
+export get_voltage_angle_series
+export get_voltage_magnitude_series
+export read_initial_conditions
+export show_states_initial_value
+export small_signal_analysis
 
 ####################################### Package Imports ####################################
 import Logging
diff --git a/src/base/jacobian.jl b/src/base/jacobian.jl
index 78a2aaa43..7b761d324 100644
--- a/src/base/jacobian.jl
+++ b/src/base/jacobian.jl
@@ -2,9 +2,12 @@
 # custom code since the structure of the system models is not compatible with the functionalities
 # in SparseDiffTools
 
-struct JacobianFunctionWrapper{F} <: Function
+struct JacobianFunctionWrapper{
+    F,
+    T <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}},
+} <: Function
     Jf::F
-    Jv::SparseArrays.SparseMatrixCSC{Float64, Int64}
+    Jv::T
     x::Vector{Float64}
 end
 
@@ -31,7 +34,7 @@ function JacobianFunctionWrapper(
     end
     Jv = SparseArrays.sparse(jac)
     Jf(Jv, x0)
-    return JacobianFunctionWrapper{typeof(Jf)}(Jf, Jv, x0)
+    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0)
 end
 
 function JacobianFunctionWrapper(
@@ -56,7 +59,7 @@ function JacobianFunctionWrapper(
     end
     Jv = SparseArrays.sparse(jac)
     Jf(Jv, x0)
-    return JacobianFunctionWrapper{typeof(Jf)}(Jf, Jv, x0)
+    return JacobianFunctionWrapper{typeof(Jf), typeof(Jv)}(Jf, Jv, x0)
 end
 
 function (J::JacobianFunctionWrapper)(x::AbstractVector{Float64})
@@ -65,33 +68,33 @@ function (J::JacobianFunctionWrapper)(x::AbstractVector{Float64})
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     x::AbstractVector{Float64},
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     J.Jf(JM, x)
     return
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     x::AbstractVector{Float64},
     p,
     t,
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     J.Jf(JM, x)
     return
 end
 
 function (J::JacobianFunctionWrapper)(
-    JM::SparseArrays.SparseMatrixCSC{Float64, Int64},
+    JM::U,
     dx::AbstractVector{Float64},
     x::AbstractVector{Float64},
     p,
     gamma,
     t,
-)
+) where {U <: Union{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}
     J.x .= x
     JM .= gamma * LinearAlgebra.Diagonal(ones(length(x))) .- J.Jf(JM, x)
     return JM
@@ -104,3 +107,23 @@ function get_jacobian(
 ) where {T <: SimulationModel}
     return JacobianFunctionWrapper(T(inputs, x0_init, JacobianCache), x0_init)
 end
+
+function _set_operating_point!(
+    x0_init::Vector{Float64},
+    inputs::SimulationInputs,
+    system::PSY.System,
+)
+    status = power_flow_solution!(x0_init, system, inputs)
+    status = initialize_static_injection!(inputs)
+    status = initialize_dynamic_injection!(x0_init, inputs, system)
+    status = initialize_dynamic_branches!(x0_init, inputs)
+    return status
+end
+function get_jacobian(T, system::PSY.System)
+    # Deepcopy avoid system modifications
+    simulation_system = deepcopy(system)
+    inputs = SimulationInputs(T, simulation_system, ReferenceBus)
+    x0_init = get_flat_start(inputs)
+    _set_operating_point!(x0_init, inputs, system)
+    return get_jacobian(T, inputs, x0_init)
+end
diff --git a/test/test_base.jl b/test/test_base.jl
index 901656b22..921fa43c0 100644
--- a/test/test_base.jl
+++ b/test/test_base.jl
@@ -504,3 +504,9 @@ end
     @test !isempty(PSY.get_components(PSY.DynamicBranch, omib_sys_copy))
     @test isempty(PSY.get_components(PSY.Line, omib_sys_copy))
 end
+
+@testset "Jacobian" begin
+    omib_sys_copy = deepcopy(omib_sys)
+    jac = get_jacobian(ResidualModel, omib_sys_copy)
+    @test isa(jac, Function)
+end

From c5c33a8ae2bb8c5af17a7f3ba054b30087e1da69 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 18:26:30 -0600
Subject: [PATCH 21/22] more test fixing

---
 test/test_case01_OMIB.jl                      | 18 ++++++-------
 test/test_case02_oneDoneQ.jl                  | 12 ++++-----
 test/test_case03_simple_marconato.jl          | 12 ++++-----
 test/test_case04_marconato.jl                 | 12 ++++-----
 test/test_case05_simple_anderson.jl           | 12 ++++-----
 test/test_case06_anderson.jl                  | 12 ++++-----
 test/test_case07_5shaft.jl                    | 12 ++++-----
 test/test_case08_VirtualSynchMachine.jl       | 12 ++++-----
 test/test_case09_oneDoneQ_inverter.jl         | 12 ++++-----
 test/test_case10_staticbranches.jl            | 12 ++++-----
 test/test_case11_dynbranches.jl               | 12 ++++-----
 test/test_case12_multimachine.jl              | 12 ++++-----
 test/test_case13_avrs.jl                      | 14 +++++-----
 test/test_case14_inverter_ref.jl              | 12 ++++-----
 test/test_case15_genrou.jl                    | 12 ++++-----
 test/test_case16_genroe.jl                    | 12 ++++-----
 test/test_case17_genrou_avr.jl                | 12 ++++-----
 test/test_case18_gensal.jl                    | 12 ++++-----
 test/test_case19_gensae.jl                    | 12 ++++-----
 test/test_case20_ac1a.jl                      | 12 ++++-----
 test/test_case21_gast.jl                      | 12 ++++-----
 test/test_case22_SteamTurbineGov1.jl          | 12 ++++-----
 test/test_case23_droopinverter.jl             | 12 ++++-----
 test/test_case24_gridfollowing.jl             | 12 ++++-----
 test/test_case25_multimach_dynlines.jl        | 12 ++++-----
 test/test_case26_SEXS.jl                      | 12 ++++-----
 test/test_case27_source_bus_voltage_change.jl | 24 ++++++++---------
 test/test_case29_renewablemodels.jl           | 26 ++++++++++---------
 test/test_case30_ieeest.jl                    | 20 ++++++++------
 test/test_case31_hygov.jl                     | 12 ++++-----
 test/x_test_sundials.jl                       |  6 ++---
 31 files changed, 207 insertions(+), 201 deletions(-)

diff --git a/test/test_case01_OMIB.jl b/test/test_case01_OMIB.jl
index 63582759d..8e74bcbca 100644
--- a/test/test_case01_OMIB.jl
+++ b/test/test_case01_OMIB.jl
@@ -32,13 +32,13 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test01_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -100,22 +100,22 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test01_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test01_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         small_sig = small_signal_analysis(sim)
         eigs = small_sig.eigenvalues
diff --git a/test/test_case02_oneDoneQ.jl b/test/test_case02_oneDoneQ.jl
index 1067d868b..5f1219035 100644
--- a/test/test_case02_oneDoneQ.jl
+++ b/test/test_case02_oneDoneQ.jl
@@ -34,13 +34,13 @@ Ybus_change = NetworkSwitch(
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test02_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -92,13 +92,13 @@ end
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test02_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case03_simple_marconato.jl b/test/test_case03_simple_marconato.jl
index 6ef944923..c652d8754 100644
--- a/test/test_case03_simple_marconato.jl
+++ b/test/test_case03_simple_marconato.jl
@@ -34,13 +34,13 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test03_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -92,13 +92,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test03_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case04_marconato.jl b/test/test_case04_marconato.jl
index 7227f9f4a..7a64db1bd 100644
--- a/test/test_case04_marconato.jl
+++ b/test/test_case04_marconato.jl
@@ -34,13 +34,13 @@ Ybus_change = NetworkSwitch(
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test04_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -91,13 +91,13 @@ end
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test04_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case05_simple_anderson.jl b/test/test_case05_simple_anderson.jl
index 76565376d..c88cd867b 100644
--- a/test/test_case05_simple_anderson.jl
+++ b/test/test_case05_simple_anderson.jl
@@ -34,13 +34,13 @@ Ybus_change = NetworkSwitch(
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test05_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -78,13 +78,13 @@ end
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test05_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case06_anderson.jl b/test/test_case06_anderson.jl
index f92aac84f..0ca594382 100644
--- a/test/test_case06_anderson.jl
+++ b/test/test_case06_anderson.jl
@@ -34,13 +34,13 @@ Ybus_change = NetworkSwitch(
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test06_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -78,13 +78,13 @@ end
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test06_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case07_5shaft.jl b/test/test_case07_5shaft.jl
index f752760cf..6d7f93464 100644
--- a/test/test_case07_5shaft.jl
+++ b/test/test_case07_5shaft.jl
@@ -37,13 +37,13 @@ Ybus_change = NetworkSwitch(
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test07_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -81,13 +81,13 @@ end
         ) #initial guess
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test07_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case08_VirtualSynchMachine.jl b/test/test_case08_VirtualSynchMachine.jl
index 110cc4b83..7cf1fedda 100644
--- a/test/test_case08_VirtualSynchMachine.jl
+++ b/test/test_case08_VirtualSynchMachine.jl
@@ -37,12 +37,12 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test08_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -93,12 +93,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test08_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case09_oneDoneQ_inverter.jl b/test/test_case09_oneDoneQ_inverter.jl
index 0c0c50866..5bffd7fff 100644
--- a/test/test_case09_oneDoneQ_inverter.jl
+++ b/test/test_case09_oneDoneQ_inverter.jl
@@ -35,12 +35,12 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))
         )
 
         # Test Initial Conditions
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -86,12 +86,12 @@ end
         )
 
         # Test Initial Conditions
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case10_staticbranches.jl b/test/test_case10_staticbranches.jl
index b5753c9e1..d02a8d5d5 100644
--- a/test/test_case10_staticbranches.jl
+++ b/test/test_case10_staticbranches.jl
@@ -36,12 +36,12 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test10_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -82,12 +82,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test10_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case11_dynbranches.jl b/test/test_case11_dynbranches.jl
index 659f7d09a..286b5289c 100644
--- a/test/test_case11_dynbranches.jl
+++ b/test/test_case11_dynbranches.jl
@@ -38,12 +38,12 @@ Ybus_change = NetworkSwitch(
         init_conds = read_initial_conditions(sim)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test10_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -78,12 +78,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test10_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case12_multimachine.jl b/test/test_case12_multimachine.jl
index 386c4ba02..814bc1c9f 100644
--- a/test/test_case12_multimachine.jl
+++ b/test/test_case12_multimachine.jl
@@ -37,13 +37,13 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test12_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -81,13 +81,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test12_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case13_avrs.jl b/test/test_case13_avrs.jl
index a79c9cc14..ec6127252 100644
--- a/test/test_case13_avrs.jl
+++ b/test/test_case13_avrs.jl
@@ -1,7 +1,7 @@
 """
 Case 13:
 This case study a three bus system with 2 machines (One d- One q-: 4th order model) and an infinite source.
-The case is similar to case 04, with different AVR and TG models.
+The case is similar to case 04, with diffvalserent AVR and TG models.
 """
 
 ##################################################
@@ -35,13 +35,13 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test13_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -76,13 +76,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test13_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case14_inverter_ref.jl b/test/test_case14_inverter_ref.jl
index 20bab6524..09ac667af 100644
--- a/test/test_case14_inverter_ref.jl
+++ b/test/test_case14_inverter_ref.jl
@@ -37,13 +37,13 @@ Ybus_change = NetworkSwitch(
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test14_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -82,13 +82,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test14_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case15_genrou.jl b/test/test_case15_genrou.jl
index d7f35f95b..3d88b412b 100644
--- a/test/test_case15_genrou.jl
+++ b/test/test_case15_genrou.jl
@@ -46,13 +46,13 @@ function test_genrou_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -106,13 +106,13 @@ function test_genrou_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case16_genroe.jl b/test/test_case16_genroe.jl
index a72550e22..495a5a4af 100644
--- a/test/test_case16_genroe.jl
+++ b/test/test_case16_genroe.jl
@@ -43,13 +43,13 @@ function test_genroe_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -103,13 +103,13 @@ function test_genroe_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case17_genrou_avr.jl b/test/test_case17_genrou_avr.jl
index 7536bad39..c41b5bae9 100644
--- a/test/test_case17_genrou_avr.jl
+++ b/test/test_case17_genrou_avr.jl
@@ -27,13 +27,13 @@ include(joinpath(TEST_FILES_DIR, "data_tests/test17.jl"))
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_genrou_avr_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -68,13 +68,13 @@ end
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_genrou_avr_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case18_gensal.jl b/test/test_case18_gensal.jl
index 2e5507f0b..5bff1e831 100644
--- a/test/test_case18_gensal.jl
+++ b/test/test_case18_gensal.jl
@@ -53,13 +53,13 @@ function test_gensal_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -113,13 +113,13 @@ function test_gensal_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case19_gensae.jl b/test/test_case19_gensae.jl
index 46a6af05c..3af6cb61a 100644
--- a/test/test_case19_gensae.jl
+++ b/test/test_case19_gensae.jl
@@ -49,13 +49,13 @@ function test_gensae_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -109,13 +109,13 @@ function test_gensae_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case20_ac1a.jl b/test/test_case20_ac1a.jl
index 635fdbd24..31d78251a 100644
--- a/test/test_case20_ac1a.jl
+++ b/test/test_case20_ac1a.jl
@@ -44,13 +44,13 @@ function test_ac1a_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -103,13 +103,13 @@ function test_ac1a_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case21_gast.jl b/test/test_case21_gast.jl
index b8a3fd86d..5c4305bad 100644
--- a/test/test_case21_gast.jl
+++ b/test/test_case21_gast.jl
@@ -27,13 +27,13 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/GAST/GAST_TEST.csv")
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_gast_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -81,13 +81,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_gast_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case22_SteamTurbineGov1.jl b/test/test_case22_SteamTurbineGov1.jl
index bb5e0003a..3bb106ab0 100644
--- a/test/test_case22_SteamTurbineGov1.jl
+++ b/test/test_case22_SteamTurbineGov1.jl
@@ -27,13 +27,13 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/TGOV1/TEST_TGOV1.csv")
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_tgov1_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -81,13 +81,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test_psse_tgov1_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case23_droopinverter.jl b/test/test_case23_droopinverter.jl
index dc16242c5..294904026 100644
--- a/test/test_case23_droopinverter.jl
+++ b/test/test_case23_droopinverter.jl
@@ -40,12 +40,12 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test23_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -93,12 +93,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test23_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case24_gridfollowing.jl b/test/test_case24_gridfollowing.jl
index 0bf47ae09..e1d2114a5 100644
--- a/test/test_case24_gridfollowing.jl
+++ b/test/test_case24_gridfollowing.jl
@@ -37,12 +37,12 @@ Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test24_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -97,12 +97,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test24_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case25_multimach_dynlines.jl b/test/test_case25_multimach_dynlines.jl
index 11682ce7d..da24bc9d5 100644
--- a/test/test_case25_multimach_dynlines.jl
+++ b/test/test_case25_multimach_dynlines.jl
@@ -36,12 +36,12 @@ Pref_change = ControlReferenceChange(1.0, gen2, :P_ref, 0.9);
         sim = Simulation!(ResidualModel, sys, path, Pref_change)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test25_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -79,12 +79,12 @@ end
         sim = Simulation!(MassMatrixModel, sys, path, Pref_change)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test25_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case26_SEXS.jl b/test/test_case26_SEXS.jl
index 8a573c57d..58f8cfeae 100644
--- a/test/test_case26_SEXS.jl
+++ b/test/test_case26_SEXS.jl
@@ -44,13 +44,13 @@ function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvalsvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvalsvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvalsvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions. Testing the simulation reset
         small_sig = small_signal_analysis(sim)
@@ -105,13 +105,13 @@ function test_sexs_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvalsvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvalsvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvalsvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions. Testing the simulation reset
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case27_source_bus_voltage_change.jl b/test/test_case27_source_bus_voltage_change.jl
index 57fd2d5ba..fa54022ce 100644
--- a/test/test_case27_source_bus_voltage_change.jl
+++ b/test/test_case27_source_bus_voltage_change.jl
@@ -37,12 +37,12 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.V_source_index,
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Solve problem
         execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02)
@@ -69,12 +69,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Solve problem
         execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02)
@@ -110,12 +110,12 @@ V_source_change = SourceBusVoltageChange(1.0, case_source, PSID.θ_source_index,
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Solve problem
         execute!(sim, IDA(), (0.0, 20.0), dtmax = 0.02)
@@ -141,12 +141,12 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test09_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Solve problem
         execute!(sim, Rodas4(), (0.0, 20.0), dtmax = 0.02)
diff --git a/test/test_case29_renewablemodels.jl b/test/test_case29_renewablemodels.jl
index 9ca21bf95..04fa6a488 100644
--- a/test/test_case29_renewablemodels.jl
+++ b/test/test_case29_renewablemodels.jl
@@ -47,13 +47,13 @@ function test_renA_implicit(csv_file, init_cond, eigs_value, F_Flag)
         sim = Simulation(ResidualModel, sys, path, Ybus_change)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -98,6 +98,7 @@ function test_renA_implicit(csv_file, init_cond, eigs_value, F_Flag)
     finally
         @info("removing test files")
         rm(path, force = true, recursive = true)
+        sim = nothing
     end
 end
 
@@ -118,13 +119,13 @@ function test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_Flag)
         sim = Simulation(MassMatrixModel, sys, path, Ybus_change)
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -169,6 +170,7 @@ function test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_Flag)
     finally
         @info("removing test files")
         rm(path, force = true, recursive = true)
+        sim = nothing
     end
 end
 
@@ -257,13 +259,13 @@ function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -319,13 +321,13 @@ function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tsp
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/test_case30_ieeest.jl b/test/test_case30_ieeest.jl
index 193cee311..aa4eafe99 100644
--- a/test/test_case30_ieeest.jl
+++ b/test/test_case30_ieeest.jl
@@ -32,7 +32,7 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
         sys = System(raw_file_dir, dyr_file)
 
         # Define Simulation Problem
-        sim = Simulation!(
+        sim = Simulation(
             ResidualModel,
             sys, #system
             path,
@@ -40,13 +40,13 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -88,7 +88,9 @@ function test_ieeest_implicit(dyr_file, csv_file, init_cond, eigs_value)
     finally
         @info("removing test files")
         rm(path, force = true, recursive = true)
+        sim = nothing
     end
+    return
 end
 
 function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
@@ -98,7 +100,7 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         sys = System(raw_file_dir, dyr_file)
 
         # Define Simulation Problem
-        sim = Simulation!(
+        sim = Simulation(
             MassMatrixModel,
             sys, #system
             path,
@@ -106,13 +108,13 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -154,7 +156,9 @@ function test_ieeest_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
     finally
         @info("removing test files")
         rm(path, force = true, recursive = true)
+        sim = nothing
     end
+    return
 end
 
 @testset "Test 30 IEEEST ResidualModel" begin
diff --git a/test/test_case31_hygov.jl b/test/test_case31_hygov.jl
index 8b616a210..be24e3ba9 100644
--- a/test/test_case31_hygov.jl
+++ b/test/test_case31_hygov.jl
@@ -27,13 +27,13 @@ csv_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/HYGOV/HYGOV_RESULTS.csv")
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test31_hygov_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
@@ -87,13 +87,13 @@ end
         )
 
         # Test Initial Condition
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test31_hygov_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions
         small_sig = small_signal_analysis(sim)
diff --git a/test/x_test_sundials.jl b/test/x_test_sundials.jl
index 998cd64ae..1ef56475e 100644
--- a/test/x_test_sundials.jl
+++ b/test/x_test_sundials.jl
@@ -55,12 +55,12 @@ function test_sundials(solver)
 
         #Obtain data for voltages
         series = get_voltage_magnitude_series(results, 102)
-        diff = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in test10_x0_init
-            diff[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
-        @test (diff[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
         @test res.solution.retcode == :Success
     finally
         @info("removing test files")

From b4ba47ed1b72930986f30b7580e3b74ac827b513 Mon Sep 17 00:00:00 2001
From: Jose Daniel Lara <jdlara@berkeley.edu>
Date: Thu, 20 Jan 2022 18:29:06 -0600
Subject: [PATCH 22/22] varname fix

---
 test/test_case26_SEXS.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/test/test_case26_SEXS.jl b/test/test_case26_SEXS.jl
index 58f8cfeae..3d73b5447 100644
--- a/test/test_case26_SEXS.jl
+++ b/test/test_case26_SEXS.jl
@@ -44,13 +44,13 @@ function test_sexs_implicit(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diffvalsvals = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diffvalsvals[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diffvalsvals[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions. Testing the simulation reset
         small_sig = small_signal_analysis(sim)
@@ -105,13 +105,13 @@ function test_sexs_mass_matrix(dyr_file, csv_file, init_cond, eigs_value)
         ) #Type of Fault
 
         # Test Initial Condition
-        diffvalsvals = [0.0]
+        diffvals = [0.0]
         res = get_init_values_for_comparison(sim)
         for (k, v) in init_cond
-            diffvalsvals[1] += LinearAlgebra.norm(res[k] - v)
+            diffvals[1] += LinearAlgebra.norm(res[k] - v)
         end
 
-        @test (diffvalsvals[1] < 1e-3)
+        @test (diffvals[1] < 1e-3)
 
         # Obtain small signal results for initial conditions. Testing the simulation reset
         small_sig = small_signal_analysis(sim)