Skip to content

Commit 325cef9

Browse files
oscardssmithoscardssmith
committed
fixed convergance
1 parent 85e5d2e commit 325cef9

File tree

2 files changed

+56
-50
lines changed

2 files changed

+56
-50
lines changed

‎lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1243,21 +1243,23 @@ end
12431243
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
12441244

12451245
# Compute linsolve_tmp for current stage
1246-
linsolve_tmp = @.. du + dtd[stage] * dT
1246+
linsolve_tmp1 = zero(du)
12471247
if mass_matrix === I
12481248
for i in 1:stage-1
1249-
linsolve_tmp = @.. linsolve_tmp + dtC[stage, i] * ks[i]
1249+
linsolve_tmp1 = @.. linsolve_tmp1 + dtC[stage, i] * ks[i]
12501250
end
12511251
else
12521252
for i in 1:stage-1
1253-
linsolve_tmp1 = mass_matrix * @..dtC[stage, i] * ks[i]
1254-
linsolve_tmp = @.. linsolve_tmp + linsolve_tmp1
1253+
linsolve_tmp1 = @.. linsolve_tmp1 + dtC[stage, i] * ks[i]
12551254
end
1255+
linsolve_tmp1 = mass_matrix * linsolve_tmp1
12561256
end
1257+
linsolve_tmp = @.. du + dtd[stage] * dT + linsolve_tmp1
12571258

12581259
ks[stage] = _reshape(W \ -_vec(linsolve_tmp), axes(uprev))
12591260
integrator.stats.nsolve += 1
12601261
end
1262+
#@show ks
12611263
u = u .+ ks[end]
12621264

12631265
if integrator.opts.adaptive
@@ -1339,19 +1341,19 @@ end
13391341
f(du, u, p, t + c[stage] * dt)
13401342
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
13411343

1344+
du1 .= 0
13421345
if mass_matrix === I
1343-
@.. linsolve_tmp = du + dtd[stage] * dT
13441346
for i in 1:stage-1
1345-
@.. linsolve_tmp += dtC[stage, i] * ks[i]
1347+
@.. du1 += dtC[stage, i] * ks[i]
13461348
end
13471349
else
1348-
du1 .= 0
13491350
for i in 1:stage-1
13501351
@.. du1 += dtC[stage, i] * ks[i]
13511352
end
13521353
mul!(_vec(du2), mass_matrix, _vec(du1))
1353-
@.. linsolve_tmp = du + dtd[stage] * dT + du2
1354+
du1 .= du2
13541355
end
1356+
@.. linsolve_tmp = du + dtd[stage] * dT + du1
13551357

13561358
linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
13571359
@.. $(_vec(ks[stage]))=-linres.u

‎lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_tableaus.jl

Lines changed: 46 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -239,32 +239,33 @@ function Rodas4Tableau(T, T2)
239239
#BET4P=0.3438D0
240240
# Coefficients for the "a" matrix
241241
A = T[
242-
0 0 0 0 0;
243-
1.544000000000000 0 0 0 0;
244-
0.9466785280815826 0.2557011698983284 0 0 0;
245-
3.314825187068521 2.896124015972201 0.9986419139977817 0 0;
246-
1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 0
242+
0 0 0 0 0 0
243+
1.544000000000000 0 0 0 0 0
244+
0.9466785280815826 0.2557011698983284 0 0 0 0
245+
3.314825187068521 2.896124015972201 0.9986419139977817 0 0 0
246+
1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 0 0
247+
1.221224509226641 6.019134481288629 12.53708332932087 -0.6878860361058950 1 0
247248
]
248249

249250
# Coefficients for the "C" matrix
250251
C = T[
251-
0 0 0 0 0;
252-
-5.668800000000000 0 0 0 0;
253-
-2.430093356833875 -0.2063599157091915 0 0 0;
254-
-0.1073529058151375 -9.594562251023355 -20.47028614809616 0 0;
255-
7.496443313967647 -10.24680431464352 -33.99990352819905 11.70890893206160 0;
252+
0 0 0 0 0
253+
-5.668800000000000 0 0 0 0
254+
-2.430093356833875 -0.2063599157091915 0 0 0
255+
-0.1073529058151375 -9.594562251023355 -20.47028614809616 0 0
256+
7.496443313967647 -10.24680431464352 -33.99990352819905 11.70890893206160 0
256257
8.083246795921522 -7.981132988064893 -31.52159432874371 16.31930543123136 -6.058818238834054
257258
]
258259

259260
# Coefficients for the "c" vector
260-
c = T2[0, 0.386, 0.21, 0.63, 1, 1]
261+
c = T2[0, 0.386, 0.21, 0.63, 1, 1, 1]
261262

262263
# Coefficients for the "d" vector
263-
d = T[0.2500000000000000, -0.1043000000000000, 0.1035000000000000, -0.03620000000000023, 0, 0]
264+
d = T[0.2500000000000000, -0.1043000000000000, 0.1035000000000000, -0.03620000000000023, 0, 0, 0]
264265

265266
# Coefficients for the "h" (kept as is)
266-
H = T[10.12623508344586 -7.487995877610167 -34.80091861555747 -7.992771707568823 1.025137723295662
267-
-0.6762803392801253 6.087714651680015 16.43084320892478 24.76722511418386 -6.594389125716872]
267+
H = T[10.12623508344586 -7.487995877610167 -34.80091861555747 -7.992771707568823 1.025137723295662 0
268+
-0.6762803392801253 6.087714651680015 16.43084320892478 24.76722511418386 -6.594389125716872 0]
268269

269270
RodasTableau(A, C, gamma, c, d, H)
270271
end
@@ -274,15 +275,16 @@ function Rodas42Tableau(T, T2)
274275

275276
# Create matrices for A and C coefficients
276277
A = T[
277-
0.0 0.0 0.0 0.0 0
278-
1.4028884 0.0 0.0 0.0 0
279-
0.6581213 -1.320936088384301 0.0 0.0 0
280-
7.131197445744498 16.02964143958207 -5.561572550509766 0.0 0
281-
22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 0
278+
0.0 0.0 0.0 0.0 0 0
279+
1.4028884 0.0 0.0 0.0 0 0
280+
0.6581213 -1.320936088384301 0.0 0.0 0 0
281+
7.131197445744498 16.02964143958207 -5.561572550509766 0.0 0 0
282+
22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 0 0
283+
22.73885722420363 67.38147284535289 -31.21877493038560 0.7285641833203814 1 0
282284
]
283285

284286
C = T[
285-
0.0 0.0 0.0 0.0 0.0 # Placeholder for C11, C12, C13, C14, C15 (unused)
287+
0.0 0.0 0.0 0.0 0.0
286288
-5.1043536 0.0 0.0 0.0 0.0
287289
-2.899967805418783 4.040399359702244 0.0 0.0 0.0
288290
-32.64449927841361 -99.35311008728094 49.99119122405989 0.0 0.0
@@ -291,12 +293,12 @@ function Rodas42Tableau(T, T2)
291293
]
292294

293295
# Create vectors for c and d coefficients
294-
c = T2[0, 0.3507221, 0.2557041, 0.6817790, 1, 1]
295-
d = T[0.2500000000000000, -0.06902209999999998, -0.0009671999999999459, -0.08797900000000025, 0,0]
296+
c = T2[0, 0.3507221, 0.2557041, 0.6817790, 1, 1, 1]
297+
d = T[0.2500000000000000, -0.06902209999999998, -0.0009671999999999459, -0.08797900000000025, 0, 0, 0]
296298

297299
# Coefficients for the "H" matrix
298-
H = T[-38.71940424117216 -135.8025833007622 64.51068857505875 -4.192663174613162 -2.531932050335060
299-
-14.99268484949843 -76.30242396627033 58.65928432851416 16.61359034616402 -0.6758691794084156]
300+
H = T[-38.71940424117216 -135.8025833007622 64.51068857505875 -4.192663174613162 -2.531932050335060 0
301+
-14.99268484949843 -76.30242396627033 58.65928432851416 16.61359034616402 -0.6758691794084156 0]
300302

301303
RodasTableau(A, C, gamma, c, d, H)
302304
end
@@ -309,11 +311,12 @@ function Rodas4PTableau(T, T2)
309311

310312
# Coefficients for the "A" matrix
311313
A = T[
312-
0.0 0.0 0.0 0.0 0
313-
3.0 0.0 0.0 0.0 0
314-
1.831036793486759 0.4955183967433795 0.0 0.0 0
315-
2.304376582692669 -0.05249275245743001 -1.176798761832782 0.0 0
316-
-7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 0
314+
0.0 0 0.0 0.0 0 0
315+
3.0 0 0.0 0.0 0 0
316+
1.831036793486759 0.4955183967433795 0.0 0.0 0 0
317+
2.304376582692669 -0.05249275245743001 -1.176798761832782 0.0 0 0
318+
-7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 0 0
319+
-7.170454962423024 -4.741636671481785 -16.31002631330971 -1.062004044111401 1 0
317320
]
318321

319322
# Coefficients for the "C" matrix
@@ -327,14 +330,14 @@ function Rodas4PTableau(T, T2)
327330
]
328331

329332
# Coefficients for the "c" vector
330-
c = T2[0, 3 * gamma, 0.21, 0.63, 1, 1]
333+
c = T2[0, 3 * gamma, 0.21, 0.63, 1, 1, 1]
331334

332335
# Coefficients for the "d" vector
333-
d = T[0.2500000000000000, -0.5000000000000000, -0.0235040000000000, -0.0362000000000000, 0, 0]
336+
d = T[0.2500000000000000, -0.5000000000000000, -0.0235040000000000, -0.0362000000000000, 0, 0, 0]
334337

335338
# Coefficients for the "H" matrix
336-
H = T[25.09876703708589 11.62013104361867 28.49148307714626 -5.664021568594133 0
337-
1.638054557396973 -0.7373619806678748 8.477918219238990 15.99253148779520 -1.882352941176471]
339+
H = T[25.09876703708589 11.62013104361867 28.49148307714626 -5.664021568594133 0 0
340+
1.638054557396973 -0.7373619806678748 8.477918219238990 15.99253148779520 -1.882352941176471 0]
338341

339342
RodasTableau(A, C, gamma, c, d, H)
340343
end
@@ -344,11 +347,12 @@ function Rodas4P2Tableau(T, T2)
344347

345348
# Coefficients for the "A" matrix
346349
A = T[
347-
0.0 0.0 0.0 0.0 0
348-
3.0 0.0 0.0 0.0 0
349-
0.906377755268814 -0.189707390391685 0.0 0.0 0
350-
3.758617027739064 1.161741776019525 -0.849258085312803 0.0 0
351-
7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 0
350+
0.0 0.0 0.0 0 0 0
351+
3.0 0.0 0.0 0 0 0
352+
0.906377755268814 -0.189707390391685 0.0 0 0 0
353+
3.758617027739064 1.161741776019525 -0.849258085312803 0 0 0
354+
7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 0 0
355+
7.089566927282776 4.573591406461604 -8.423496976860259 -0.959280113459775 1 0
352356
]
353357

354358
# Coefficients for the "C" matrix
@@ -362,14 +366,14 @@ function Rodas4P2Tableau(T, T2)
362366
]
363367

364368
# Coefficients for the "c" vector
365-
c = T2[0, 0.750000000000000, 0.321448134013046, 0.519745732277726, 1, 1]
369+
c = T2[0, 0.750000000000000, 0.321448134013046, 0.519745732277726, 1, 1, 1]
366370

367371
# Coefficients for the "d" vector
368-
d = T[0.250000000000000, -0.500000000000000, -0.189532918363016, 0.085612108792769, 0, 0]
372+
d = T[0.250000000000000, -0.500000000000000, -0.189532918363016, 0.085612108792769, 0, 0, 0]
369373

370374
# Coefficients for the "H" matrix
371-
H = [-5.323528268423303 -10.042123754867493 17.175254928256965 -5.079931171878093 -0.016185991706112
372-
6.984505741529879 6.914061169603662 -0.849178943070653 18.104410789349338 -3.516963011559032]
375+
H = [-5.323528268423303 -10.042123754867493 17.175254928256965 -5.079931171878093 -0.016185991706112 0
376+
6.984505741529879 6.914061169603662 -0.849178943070653 18.104410789349338 -3.516963011559032 0]
373377

374378
RodasTableau(A, C, gamma, c, d, H)
375379
end

0 commit comments

Comments
 (0)