Skip to content

Commit df027cf

Browse files
claudeclaude
committed
Complete NAND Gate DAE benchmark following TransistorAmplifier pattern
- Add DAEProblem formulation using ModelingToolkit conversion - Add ModelingToolkit version via modelingtoolkitize and structural_simplify - Expand solver testing to include CVODE_BDF, DFBDF, DASKR, radau, RadauIIA5 - Implement three problem formulations: MTK, DAE, and mass matrix versions - Add comprehensive work-precision benchmarks at multiple tolerance levels - Include timeseries error analysis for all formulations - Follow established benchmark structure with prob_choice parameter - Add required imports: Sundials, DASSL, DASKR 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <[email protected]>
1 parent 2108701 commit df027cf

File tree

1 file changed

+146
-91
lines changed

1 file changed

+146
-91
lines changed

‎benchmarks/DAE/NANDGateProblem.jmd‎

Lines changed: 146 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ author: Jayant Pranjal
55

66
```julia
77
using OrdinaryDiffEq, DiffEqDevTools, ModelingToolkit, ODEInterfaceDiffEq,
8-
Plots
8+
Plots, Sundials, DASSL, DASKR
99
using LinearAlgebra
1010
using ModelingToolkit: t_nounits as t, D_nounits as D
1111
```
@@ -151,31 +151,6 @@ function ibd(vbd)
151151
end
152152
```
153153

154-
## Capacitance Matrix Function
155-
156-
```julia
157-
function capacitance_matrix!(C, y, t)
158-
fill!(C, 0.0)
159-
160-
C[1,1] = CGS
161-
C[2,2] = CGD
162-
C[3,3] = CBS
163-
C[4,4] = CBD
164-
C[5,5] = 0.0
165-
C[6,6] = CGS
166-
C[7,7] = CGD
167-
C[8,8] = CBS
168-
C[9,9] = CBD
169-
C[10,10] = 0.0
170-
C[11,11] = CGS
171-
C[12,12] = CGD
172-
C[13,13] = CBS
173-
C[14,14] = CBD
174-
175-
return C
176-
end
177-
```
178-
179154
## DAE System Definition
180155

181156
```julia
@@ -207,121 +182,201 @@ function nand_rhs!(f, y, p, t)
207182
return nothing
208183
end
209184

210-
function nand_mass_matrix!(M, y, p, t)
211-
capacitance_matrix!(M, y, t)
212-
return nothing
213-
end
185+
# Capacitance matrix (mass matrix)
186+
dirMassMatrix = [
187+
CGS 0 0 0 0 0 0 0 0 0 0 0 0 0
188+
0 CGD 0 0 0 0 0 0 0 0 0 0 0 0
189+
0 0 CBS 0 0 0 0 0 0 0 0 0 0 0
190+
0 0 0 CBD 0 0 0 0 0 0 0 0 0 0
191+
0 0 0 0 0 0 0 0 0 0 0 0 0 0
192+
0 0 0 0 0 CGS 0 0 0 0 0 0 0 0
193+
0 0 0 0 0 0 CGD 0 0 0 0 0 0 0
194+
0 0 0 0 0 0 0 CBS 0 0 0 0 0 0
195+
0 0 0 0 0 0 0 0 CBD 0 0 0 0 0
196+
0 0 0 0 0 0 0 0 0 0 0 0 0 0
197+
0 0 0 0 0 0 0 0 0 0 CGS 0 0 0
198+
0 0 0 0 0 0 0 0 0 0 0 CGD 0 0
199+
0 0 0 0 0 0 0 0 0 0 0 0 CBS 0
200+
0 0 0 0 0 0 0 0 0 0 0 0 0 CBD
201+
]
214202

215-
function create_nand_problem()
216-
y0 = [5.0, 5.0, VBB, VBB, 5.0, 3.62385, 5.0, VBB, VBB, 3.62385, 0.0, 3.62385, VBB, VBB]
217-
dy0 = zeros(14)
218-
tspan = (0.0, 80.0)
219-
M = zeros(14, 14)
220-
capacitance_matrix!(M, y0, 0.0)
221-
f = ODEFunction(nand_rhs!, mass_matrix=M)
222-
prob = ODEProblem(f, y0, tspan)
223-
return prob
224-
end
225-
nand_prob = create_nand_problem()
203+
# Initial conditions
204+
y0 = [5.0, 5.0, VBB, VBB, 5.0, 3.62385, 5.0, VBB, VBB, 3.62385, 0.0, 3.62385, VBB, VBB]
205+
tspan = (0.0, 80.0)
206+
207+
# Create the three problem formulations
208+
mmf = ODEFunction(nand_rhs!, mass_matrix=dirMassMatrix)
209+
mmprob = ODEProblem(mmf, y0, tspan)
210+
211+
# ModelingToolkit version via modelingtoolkitize
212+
mtk_sys = modelingtoolkitize(mmprob)
213+
mtk_simplified = structural_simplify(mtk_sys)
214+
u0_mtk = [mtk_simplified.states[i] => mmprob.u0[i] for i in 1:length(mmprob.u0)]
215+
mtkprob = ODEProblem(mtk_simplified, u0_mtk, mmprob.tspan)
216+
217+
# DAEProblem version
218+
du = mmprob.f(mmprob.u0, mmprob.p, 0.0)
219+
du0 = D.(unknowns(mtk_simplified)) .=> du
220+
daeprob = DAEProblem(mtk_simplified, du0, [], tspan)
221+
222+
# Generate reference solutions
223+
ref_sol = solve(mtkprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
224+
dae_ref_sol = solve(daeprob, IDA(), abstol=1e-10, reltol=1e-10, tstops=0.0:5.0:80.0)
225+
mm_ref_sol = solve(mmprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
226+
227+
probs = [mtkprob, daeprob, mmprob]
228+
refs = [ref_sol, ref_sol, mm_ref_sol]
226229
```
227230

228-
## Generate Reference Solution
231+
## Generate Reference Solution and Plot
229232

230233
```julia
231-
ref_sol = solve(nand_prob, Rodas5P(), abstol=1e-12, reltol=1e-12,
232-
tstops=0.0:5.0:80.0) # Include discontinuity points
234+
plot(ref_sol, title="NAND Gate Circuit - Node Potentials (MTK)",
235+
xlabel="Time", ylabel="Voltage (V)", legend=:outertopright)
236+
```
233237

234-
plot(ref_sol, title="NAND Gate Circuit - Node Potentials",
238+
```julia
239+
plot(mm_ref_sol, title="NAND Gate Circuit - Node Potentials (Mass Matrix)",
235240
xlabel="Time", ylabel="Voltage (V)", legend=:outertopright)
236241
```
237242

243+
## Omissions
244+
245+
Some solvers are omitted due to performance or stability issues at certain tolerances.
246+
238247
## Work-Precision Benchmarks
239248

240249
### High Tolerances
241250

242251
```julia
243-
abstols = 1.0 ./ 10.0 .^ (4:7)
252+
abstols = 1.0 ./ 10.0 .^ (5:8)
244253
reltols = 1.0 ./ 10.0 .^ (1:4)
245254

246255
setups = [
247-
Dict(:alg=>Rosenbrock23()),
248-
Dict(:alg=>Rodas4()),
249-
Dict(:alg=>Rodas5P()),
250-
Dict(:alg=>FBDF()),
251-
Dict(:alg=>QNDF()),
252-
253-
256+
Dict(:prob_choice => 1, :alg=>Rodas4()),
257+
Dict(:prob_choice => 1, :alg=>FBDF()),
258+
Dict(:prob_choice => 1, :alg=>QNDF()),
259+
Dict(:prob_choice => 1, :alg=>radau()),
260+
Dict(:prob_choice => 1, :alg=>RadauIIA5()),
261+
Dict(:prob_choice => 2, :alg=>DFBDF()),
262+
Dict(:prob_choice => 2, :alg=>IDA()),
263+
Dict(:prob_choice => 2, :alg=>CVODE_BDF())
254264
]
255265

256-
wp = WorkPrecisionSet([nand_prob], abstols, reltols, setups;
257-
save_everystep=false, appxsol=[ref_sol],
258-
maxiters=Int(1e6), numruns=5,
266+
wp = WorkPrecisionSet(probs, abstols, reltols, setups;
267+
save_everystep=false, appxsol=refs,
268+
maxiters=Int(1e5), numruns=10,
259269
tstops=0.0:5.0:80.0)
260270
plot(wp, title="NAND Gate DAE - Work-Precision (High Tolerances)")
261271
```
262272

263-
### Medium Tolerances
264-
265273
```julia
266-
abstols = 1.0 ./ 10.0 .^ (6:9)
267-
reltols = 1.0 ./ 10.0 .^ (3:6)
274+
abstols = 1.0 ./ 10.0 .^ (6:8)
275+
reltols = 1.0 ./ 10.0 .^ (2:4)
268276

269277
setups = [
270-
Dict(:alg=>Rodas4()),
271-
Dict(:alg=>Rodas5P()),
272-
Dict(:alg=>FBDF()),
273-
Dict(:alg=>QNDF()),
274-
275-
278+
Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
279+
Dict(:prob_choice => 1, :alg=>Rodas4()),
280+
Dict(:prob_choice => 2, :alg=>IDA()),
281+
Dict(:prob_choice => 3, :alg=>Rodas5P()),
282+
Dict(:prob_choice => 3, :alg=>Rodas4()),
283+
Dict(:prob_choice => 3, :alg=>FBDF()),
284+
Dict(:prob_choice => 2, :alg=>IDA()),
285+
Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
286+
Dict(:prob_choice => 2, :alg=>DASKR.daskr())
276287
]
277288

278-
wp = WorkPrecisionSet([nand_prob], abstols, reltols, setups;
279-
save_everystep=false, appxsol=[ref_sol],
280-
maxiters=Int(1e6), numruns=5,
289+
wp = WorkPrecisionSet(probs, abstols, reltols, setups;
290+
save_everystep=false, appxsol=refs,
291+
maxiters=Int(1e5), numruns=10,
281292
tstops=0.0:5.0:80.0)
282293
plot(wp, title="NAND Gate DAE - Work-Precision (Medium Tolerances)")
283294
```
284295

285-
### Low Tolerances (High Accuracy)
296+
### Timeseries Errors
286297

287298
```julia
288-
abstols = 1.0 ./ 10.0 .^ (8:11)
289-
reltols = 1.0 ./ 10.0 .^ (5:8)
299+
abstols = 1.0 ./ 10.0 .^ (5:8)
300+
reltols = 1.0 ./ 10.0 .^ (1:4)
290301

291302
setups = [
292-
Dict(:alg=>Rodas5P()),
293-
Dict(:alg=>FBDF()),
294-
Dict(:alg=>QNDF()),
303+
Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
304+
Dict(:prob_choice => 1, :alg=>Rodas4()),
305+
Dict(:prob_choice => 1, :alg=>FBDF()),
306+
Dict(:prob_choice => 1, :alg=>QNDF()),
307+
Dict(:prob_choice => 1, :alg=>radau()),
308+
Dict(:prob_choice => 1, :alg=>RadauIIA5()),
309+
Dict(:prob_choice => 2, :alg=>DFBDF()),
310+
Dict(:prob_choice => 2, :alg=>IDA()),
311+
Dict(:prob_choice => 2, :alg=>CVODE_BDF())
312+
]
313+
314+
wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate=:l2,
315+
save_everystep=false, appxsol=refs,
316+
maxiters=Int(1e5), numruns=10,
317+
tstops=0.0:5.0:80.0)
318+
plot(wp, title="NAND Gate DAE - Timeseries Errors (High Tolerances)")
319+
```
295320

321+
```julia
322+
abstols = 1.0 ./ 10.0 .^ (6:8)
323+
reltols = 1.0 ./ 10.0 .^ (2:4)
296324

325+
setups = [
326+
Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
327+
Dict(:prob_choice => 1, :alg=>Rodas4()),
328+
Dict(:prob_choice => 2, :alg=>IDA()),
329+
Dict(:prob_choice => 3, :alg=>Rodas5P()),
330+
Dict(:prob_choice => 3, :alg=>Rodas4()),
331+
Dict(:prob_choice => 3, :alg=>FBDF()),
332+
Dict(:prob_choice => 2, :alg=>IDA()),
333+
Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
334+
Dict(:prob_choice => 2, :alg=>DASKR.daskr())
297335
]
298336

299-
wp = WorkPrecisionSet([nand_prob], abstols, reltols, setups;
300-
save_everystep=false, appxsol=[ref_sol],
301-
maxiters=Int(1e6), numruns=5,
337+
wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate=:l2,
338+
save_everystep=false, appxsol=refs,
339+
maxiters=Int(1e5), numruns=10,
302340
tstops=0.0:5.0:80.0)
303-
plot(wp, title="NAND Gate DAE - Work-Precision (Low Tolerances)")
341+
plot(wp, title="NAND Gate DAE - Timeseries Errors (Medium Tolerances)")
304342
```
305343

306-
### Timeseries Error Analysis
344+
### Low Tolerances (High Accuracy)
345+
346+
This measures performance when high accuracy is needed.
307347

308348
```julia
309-
abstols = 1.0 ./ 10.0 .^ (5:8)
310-
reltols = 1.0 ./ 10.0 .^ (2:5)
349+
abstols = 1.0 ./ 10.0 .^ (7:12)
350+
reltols = 1.0 ./ 10.0 .^ (4:9)
311351

312352
setups = [
313-
Dict(:alg=>Rodas4()),
314-
Dict(:alg=>Rodas5P()),
315-
Dict(:alg=>FBDF()),
316-
317-
353+
Dict(:prob_choice => 1, :alg=>Rodas5P()),
354+
Dict(:prob_choice => 3, :alg=>Rodas5P()),
355+
Dict(:prob_choice => 1, :alg=>Rodas4()),
356+
Dict(:prob_choice => 3, :alg=>Rodas4()),
357+
Dict(:prob_choice => 1, :alg=>FBDF()),
358+
Dict(:prob_choice => 1, :alg=>QNDF()),
359+
Dict(:prob_choice => 1, :alg=>radau()),
360+
Dict(:prob_choice => 1, :alg=>RadauIIA5()),
361+
Dict(:prob_choice => 2, :alg=>DFBDF()),
362+
Dict(:prob_choice => 2, :alg=>IDA()),
363+
Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
364+
Dict(:prob_choice => 2, :alg=>DASKR.daskr())
318365
]
319366

320-
wp = WorkPrecisionSet([nand_prob], abstols, reltols, setups;
321-
error_estimate=:l2, save_everystep=false,
322-
appxsol=[ref_sol], maxiters=Int(1e6), numruns=5,
367+
wp = WorkPrecisionSet(probs, abstols, reltols, setups;
368+
save_everystep=false, appxsol=refs,
369+
maxiters=Int(1e5), numruns=10,
370+
tstops=0.0:5.0:80.0)
371+
plot(wp, title="NAND Gate DAE - Work-Precision (Low Tolerances)")
372+
```
373+
374+
```julia
375+
wp = WorkPrecisionSet(probs, abstols, reltols, setups; error_estimate=:l2,
376+
save_everystep=false, appxsol=refs,
377+
maxiters=Int(1e5), numruns=10,
323378
tstops=0.0:5.0:80.0)
324-
plot(wp, title="NAND Gate DAE - Timeseries Error Analysis")
379+
plot(wp, title="NAND Gate DAE - Timeseries Errors (Low Tolerances)")
325380
```
326381

327382
## Analysis of Key Nodes

0 commit comments

Comments
 (0)