remove ODEProblem and associated structure

src/HallThruster.jl

@@ -54,7 +54,6 @@ include("collisions/collision_frequencies.jl")
 
 include("simulation/initialization.jl")
 include("simulation/geometry.jl")
-include("simulation/postprocess.jl")
 include("simulation/boundaryconditions.jl")
 include("simulation/potential.jl")
 include("simulation/heavy_species.jl")
@@ -63,9 +62,10 @@ include("simulation/sourceterms.jl")
 include("simulation/plume.jl")
 include("simulation/update_electrons.jl")
 include("simulation/configuration.jl")
-include("simulation/restart.jl")
-include("simulation/ssprk_step.jl")
+include("simulation/solution.jl")
 include("simulation/simulation.jl")
+include("simulation/restart.jl")
+include("simulation/postprocess.jl")
 include("visualization/plotting.jl")
 include("visualization/recipes.jl")
 

src/simulation/postprocess.jl

@@ -1,21 +1,3 @@
-struct Solution{T, U, P, S}
-    t::T
-    u::U
-    savevals::S
-    retcode::Symbol
-    params::P
-end
-
-function Solution(sol::S, params::P, savevals::SV) where {S, P, SV}
-    return Solution(sol.t, sol.u, savevals, sol.retcode, params)
-end
-
-function Base.show(io::IO, mime::MIME"text/plain", sol::Solution)
-    println(io, "Hall thruster solution with $(length(sol.u)) saved frames")
-    println(io, "Retcode: $(string(sol.retcode))")
-    print(io, "End time: $(sol.t[end]) seconds")
-end
-
 """
     time_average(sol, tstampstart)
 compute time-averaged solution, input Solution type and the frame at which averaging starts.

src/simulation/simulation.jl

@@ -249,9 +249,8 @@ function run_simulation(
     # make values in params available for first timestep
     update_electrons!(U, params)
 
-    # Set up and solve problem, this form is temporary
-    prob = ODEProblem(update_heavy_species!, U, tspan, params)
-    sol_info = @timed solve(prob; saveat)
+    # Set up
+    sol_info = @timed solve(U, params, tspan; saveat)
 
     sol = sol_info.value
     sim_time = sol_info.time

src/simulation/ssprk_step.jl → src/simulation/solution.jl

@@ -1,30 +1,39 @@
 # Perform one step of the Strong-stability-preserving RK22 algorithm
-function ssprk22_step!(u, f, params, t, dt)
-    (;cache) = params
-    (;k, u1) = cache
+function ssprk22_step!(U, f, params, t, dt)
+    (;k, u1) = params.cache
 
     # First step of SSPRK22
-    f(k, u, params, t)
-    @. u1 = u + dt * k
+    f(k, U, params, t)
+    @. u1 = U + dt * k
     stage_limiter!(u1, params)
 
     # Second step of SSPRK22
     f(k, u1, params, t+dt)
-    @. u = (u + u1 + dt * k) / 2
-    stage_limiter!(u, params)
+    @. U = (U + u1 + dt * k) / 2
+    stage_limiter!(U, params)
 
     return nothing
 end
 
-struct ODEProblem{F, U, T, P}
-    f::F
+struct Solution{T, U, P, S}
+    t::T
     u::U
-    tspan::T
+    savevals::S
+    retcode::Symbol
     params::P
 end
 
-function solve(prob::ODEProblem; saveat)
-    (;f, u, tspan, params) = prob
+function Solution(sol::S, params::P, savevals::SV) where {S, P, SV}
+    return Solution(sol.t, sol.u, savevals, sol.retcode, params)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain", sol::Solution)
+    println(io, "Hall thruster solution with $(length(sol.u)) saved frames")
+    println(io, "Retcode: $(string(sol.retcode))")
+    print(io, "End time: $(sol.t[end]) seconds")
+end
+
+function solve(U, params, tspan; saveat)
     i = 1
     save_ind = 2
     t = tspan[1]
@@ -40,10 +49,10 @@ function solve(prob::ODEProblem; saveat)
     )
 
     first_saveval = NamedTuple{fields_to_save}(params.cache)
-    u_save = [deepcopy(u) for _ in saveat]
+    u_save = [deepcopy(U) for _ in saveat]
     savevals = [deepcopy(first_saveval) for _ in saveat]
 
-    (nvars, ncells) = size(u)
+    (nvars, ncells) = size(U)
 
     while t < tspan[2]
         i += 1
@@ -54,35 +63,36 @@ function solve(prob::ODEProblem; saveat)
         t += params.dt[]
 
         # Update heavy species
-        ssprk22_step!(u, f, params, t, params.dt[])
+        ssprk22_step!(U, update_heavy_species!, params, t, params.dt[])
 
         # Update electron quantities
-        update_electrons!(u, params, t)
+        update_electrons!(U, params, t)
 
         # Update plume geometry
-        update_plume_geometry!(u, params)
+        update_plume_geometry!(U, params)
 
         # Check for NaNs and terminate if necessary
         nandetected = false
         infdetected = false
 
         @inbounds for j in 1:ncells, i in 1:nvars
-            if isnan(u[i, j])
+            if isnan(U[i, j])
                 @warn("NaN detected in variable $i in cell $j at time $(t)")
                 nandetected = true
                 retcode = :NaNDetected
                 break
-            elseif isinf(u[i, j])
+            elseif isinf(U[i, j])
                 @warn("Inf detected in variable $i in cell $j at time $(t)")
                 infdetected = true
                 retcode = :InfDetected
                 break
             end
         end
 
-        # Save values, TODO interpolate these to be exact
+        # Save values at designated intervals
+        # TODO interpolate these to be exact
         if t > saveat[save_ind]
-            u_save[save_ind] .= u
+            u_save[save_ind] .= U
             for field in fields_to_save
                 if field == :anom_variables
                     for i in 1:num_anom_variables(params.config.anom_model)

src/utilities/utility_functions.jl

@@ -78,7 +78,6 @@ function bitwise_log2ceil(x)
 end
 
 @inline function find_left_index(val, r::LinRangeWrapper)
-    #ceil(Int, clamp(muladd(r.A, val, r.B), 0, r.r.len))
     return searchsortedfirst(r.r, val) - 1
 end
 
@@ -120,7 +119,6 @@ function cumtrapz(x, y, y0 = zero(typeof(y[1] * x[1])))
 end
 
 function cumtrapz!(cache, x, y, y0 = zero(typeof(y[1] * x[1])))
-
     cache[1] = y0
     @inbounds for i in 2:lastindex(x)
         Δx = x[i] - x[i-1]

test/runtests.jl

@@ -24,9 +24,10 @@ include("unit_tests/test_restarts.jl")
 include("unit_tests/test_json.jl")
 
 using Symbolics
+include("order_verification/ovs_funcs.jl")
+include("order_verification/ovs_energy.jl")
+include("order_verification/ovs_ions.jl")
 @testset "Order verification (electron energy)" begin
-    include("order_verification/ovs_funcs.jl")
-    include("order_verification/ovs_energy.jl")
 
     vfunc = x -> OVS_Energy.verify_energy(x)
     refinements = refines(6, 20, 2)
@@ -48,8 +49,6 @@ using Symbolics
 end
 
 @testset "Order verification (neutrals and ions)" begin
-    include("order_verification/ovs_funcs.jl")
-    include("order_verification/ovs_ions.jl")
     refinements = refines(6, 10, 2)
 
     limiter = HallThruster.van_leer