@@ -208,25 +208,51 @@ tspan = (0.0, 80.0)
208208mmf = ODEFunction(nand_rhs!, mass_matrix=dirMassMatrix)
209209mmprob = ODEProblem(mmf, y0, tspan)
210210
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 = zeros(length(mmprob.u0))
219- mmprob.f(du, mmprob.u0, mmprob.p, 0.0)
220- du0 = D.(unknowns(mtk_simplified)) .=> du
221- daeprob = DAEProblem(mtk_simplified, du0, [], tspan)
211+ # DAEProblem version using direct DAE formulation
212+ function nand_dae!(out, du, u, p, t)
213+ v1 = V1(t)
214+ v2 = V2(t)
215+ v1d = V1_derivative(t)
216+ v2d = V2_derivative(t)
217+
218+ y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14 = u
219+ dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8, dy9, dy10, dy11, dy12, dy13, dy14 = du
220+
221+ # Differential equations (with capacitors)
222+ out[1] = CGS * dy1 + (y1 - y5) / RGS + ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
223+ out[2] = CGD * dy2 + (y2 - VDD) / RGD - ids(1, y2 - y1, y5 - y1, y3 - y5, y5 - y2, y4 - VDD)
224+ out[3] = CBS * dy3 + (y3 - VBB) / RBS - ibs(y3 - y5)
225+ out[4] = CBD * dy4 + (y4 - VBB) / RBD - ibd(y4 - VDD)
226+
227+ # Algebraic equations (KCL at nodes)
228+ out[5] = (y5 - y1) / RGS + ibs(y3 - y5) + (y5 - y7) / RGD + ibd(y9 - y5)
229+
230+ out[6] = CGS * dy6 - CGS * v1d + (y6 - y10) / RGS + ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
231+ out[7] = CGD * dy7 - CGD * v1d + (y7 - y5) / RGD - ids(2, y7 - y6, v1 - y6, y8 - y10, v1 - y7, y9 - y5)
232+ out[8] = CBS * dy8 + (y8 - VBB) / RBS - ibs(y8 - y10)
233+ out[9] = CBD * dy9 + (y9 - VBB) / RBD - ibd(y9 - y5)
234+
235+ # Algebraic equation (KCL at node)
236+ out[10] = (y10 - y6) / RGS + ibs(y8 - y10) + (y10 - y12) / RGD + ibd(y14 - y10)
237+
238+ out[11] = CGS * dy11 - CGS * v2d + y11 / RGS + ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
239+ out[12] = CGD * dy12 - CGD * v2d + (y12 - y10) / RGD - ids(2, y12 - y11, v2 - y11, y13, v2 - y12, y14 - y10)
240+ out[13] = CBS * dy13 + (y13 - VBB) / RBS - ibs(y13)
241+ out[14] = CBD * dy14 + (y14 - VBB) / RBD - ibd(y14 - y10)
242+
243+ return nothing
244+ end
245+
246+ # Create DAE problem
247+ du0_dae = zeros(14)
248+ daeprob = DAEProblem(nand_dae!, du0_dae, y0, tspan)
222249
223250# Generate reference solutions
224251ref_sol = solve(mmprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
225252dae_ref_sol = solve(daeprob, IDA(), abstol=1e-10, reltol=1e-10, tstops=0.0:5.0:80.0)
226- mtk_ref_sol = solve(mtkprob, Rodas5P(), abstol=1e-12, reltol=1e-12, tstops=0.0:5.0:80.0)
227253
228- probs = [mmprob, daeprob, mtkprob ]
229- refs = [ref_sol, dae_ref_sol, mtk_ref_sol ]
254+ probs = [mmprob, daeprob]
255+ refs = [ref_sol, dae_ref_sol]
230256```
231257
232258## Generate Reference Solution and Plot
@@ -260,7 +286,8 @@ setups = [
260286 Dict(:prob_choice => 1, :alg=>radau()),
261287 Dict(:prob_choice => 1, :alg=>RadauIIA5()),
262288 Dict(:prob_choice => 2, :alg=>IDA()),
263- Dict(:prob_choice => 2, :alg=>CVODE_BDF())
289+ Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
290+ Dict(:prob_choice => 2, :alg=>DASKR.daskr())
264291]
265292
266293wp = WorkPrecisionSet(probs, abstols, reltols, setups;
@@ -277,10 +304,8 @@ reltols = 1.0 ./ 10.0 .^ (2:4)
277304setups = [
278305 Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
279306 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()),
307+ Dict(:prob_choice => 1, :alg=>Rodas5P()),
308+ Dict(:prob_choice => 1, :alg=>FBDF()),
284309 Dict(:prob_choice => 2, :alg=>IDA()),
285310 Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
286311 Dict(:prob_choice => 2, :alg=>DASKR.daskr())
@@ -324,10 +349,8 @@ reltols = 1.0 ./ 10.0 .^ (2:4)
324349setups = [
325350 Dict(:prob_choice => 1, :alg=>Rosenbrock23()),
326351 Dict(:prob_choice => 1, :alg=>Rodas4()),
327- Dict(:prob_choice => 2, :alg=>IDA()),
328- Dict(:prob_choice => 3, :alg=>Rodas5P()),
329- Dict(:prob_choice => 3, :alg=>Rodas4()),
330- Dict(:prob_choice => 3, :alg=>FBDF()),
352+ Dict(:prob_choice => 1, :alg=>Rodas5P()),
353+ Dict(:prob_choice => 1, :alg=>FBDF()),
331354 Dict(:prob_choice => 2, :alg=>IDA()),
332355 Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
333356 Dict(:prob_choice => 2, :alg=>DASKR.daskr())
@@ -350,14 +373,11 @@ reltols = 1.0 ./ 10.0 .^ (4:9)
350373
351374setups = [
352375 Dict(:prob_choice => 1, :alg=>Rodas5P()),
353- Dict(:prob_choice => 3, :alg=>Rodas5P()),
354376 Dict(:prob_choice => 1, :alg=>Rodas4()),
355- Dict(:prob_choice => 3, :alg=>Rodas4()),
356377 Dict(:prob_choice => 1, :alg=>FBDF()),
357378 Dict(:prob_choice => 1, :alg=>QNDF()),
358379 Dict(:prob_choice => 1, :alg=>radau()),
359380 Dict(:prob_choice => 1, :alg=>RadauIIA5()),
360- Dict(:prob_choice => 2, :alg=>DFBDF()),
361381 Dict(:prob_choice => 2, :alg=>IDA()),
362382 Dict(:prob_choice => 2, :alg=>CVODE_BDF()),
363383 Dict(:prob_choice => 2, :alg=>DASKR.daskr())
0 commit comments