ODE Solver Benchmarking

.gitignore

@@ -9,6 +9,8 @@
 .vscode
 *.DS_Store
 *.pdf
+*.png
+*.jld2
 .gitattributes
 Manifest.toml
 build.log

benchmarks/DIIID_ideal_example/dcon.toml

@@ -5,14 +5,15 @@ ode_flag = true              # Integrate ODE's for determining stability of inte
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 mer_flag = true               # Evaluate the Mercier criterian
 
-set_psilim_via_dmlim = true               # Safety factor (q) limit determined as q_ir+dmlim...
-dmlim = 0.2                   # See sas_flag
+set_psilim_via_dmlim = true   # Safety factor (q) limit determined as q_ir+dmlim...
+dmlim = 0.2                   # See set_psilim_via_dmlim
 psiedge = 1.00                # If less then psilim, calculates dW(psi)...
 qlow = 1.02                   # Integration initiated at q determined by min(q0, qlow)...
 qhigh = 1e3                   # Integration terminated at q limit determined by min(qa, qhigh)...
 sing_start = 0                # Start integration at the sing_start'th rational from the axis (psilow)
 
-nn = 1                        # Toroidal mode number
+nn_low = 1                    # Smallest toroidal mode number to include
+nn_high = 1                   # Largest toroidal mode number to include
 delta_mlow = 8                # Expands lower bound of Fourier harmonics
 delta_mhigh = 8               # Expands upper bound of Fourier harmonics
 delta_mband = 0               # Integration keeps only this wide a band...
@@ -37,5 +38,15 @@ crossover = 1e-2              # Fractional distance from rational q at which tol
 singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
 ucrit = 1e4                   # Maximum fraction of solutions allowed before re-normalized
 
+[WALL]
+shape = "nowall"
+aw = 0.05
+bw = 1.5
+cw = 0
+dw = 0.5
+tw = 0.05
+equal_arc_wall = 0
+a = 0.2415
+
 [DEBUG]
 output_benchmark_data = true      # Output benchmark data for comparison between codes

src/DCON/Main.jl

@@ -176,6 +176,8 @@ function Main(path::String="./")
     println("Normal termination.")
 
     # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found
+
+    return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, vac_data=ctrl.vac_flag ? vac_data : nothing)
 end
 
 """

src/DCON/Ode.jl

@@ -298,7 +298,7 @@ function ode_step!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaE
     # Advance differential equation to next singular surface or edge
     rtol = compute_tols(ctrl, intr, odet) # initial tolerances
     prob = ODEProblem(sing_der!, odet.u, (odet.psifac, odet.psimax), (ctrl, equil, ffit, intr, odet))
-    sol = solve(prob, Tsit5(); reltol=rtol, callback=cb)
+    sol = solve(prob, BS5(); reltol=rtol, callback=cb)
     # TODO: check absolute tolerances, check how sensitive outputs are to tolerances
 
     # Update u and psifac with the solution at the end of the interval

src/Equilibrium/DirectEquilibrium.jl

@@ -294,7 +294,7 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false))
 
     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
-    sol = solve(prob, Tsit5(); callback=callback, reltol=1e-6, abstol=1e-8, dt=2π / 200, adaptive=true)
+    sol = solve(prob, BS5(); callback=callback, reltol=1e-6, abstol=1e-8, dt=2π / 200, adaptive=true)
 
     if sol.retcode != :Success && sol.retcode != :Terminated
         error("ODE integration failed for psi = $psifac with code: $(sol.retcode)")

test/runtests_fullruns.jl

@@ -2,9 +2,15 @@
 @testset "Full DCON runs" begin
     ex1 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
     @info "Running Solovev ideal example"
-    @test isnothing(DCON.Main(ex1))
-
+    @test begin
+        DCON.Main(ex1)
+        true
+    end
+
     ex2 = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example_multi_n")
     @info "Running Solovev ideal multi-n example"
-    @test isnothing(DCON.Main(ex2))
+    @test begin
+        DCON.Main(ex2)
+        true
+    end
 end