Skip to content

Commit 6e494cc

Browse files
ivan-piivan-pi
committed
Add Caroll's test cases, Hosea & Shampine's test cases and TRBDF2 solver
1 parent c82191d commit 6e494cc

7 files changed

Lines changed: 1123 additions & 1 deletion

File tree

‎example/CMakeLists.txt‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,12 @@ endforeach()
77
# The test takes too long
88
add_executable(ring_modulator ring_modulator.f90)
99
target_link_libraries(ring_modulator PRIVATE stiff3)
10+
11+
add_executable(caroll_driver caroll_driver.f90 caroll_test_problems.f90)
12+
target_link_libraries(caroll_driver PRIVATE stiff3)
13+
14+
add_executable(hosea_driver hosea_driver.f90)
15+
target_link_libraries(hosea_driver PRIVATE stiff3)
16+
17+
add_executable(fluid_bed_reactor fluid_bed_reactor.f90)
18+
target_link_libraries(fluid_bed_reactor PRIVATE stiff3)

‎example/caroll_driver.f90‎

Lines changed: 261 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,261 @@
1+
! caroll_driver.f90 -- A set of of ODE test problems
2+
!
3+
! This driver implements the test problems used in
4+
!
5+
! Caroll, John (1989). A composite integration scheme for the numerical
6+
! solution of systems of ordinary differential equations.
7+
! Journal of Computational and Applied Mathematics 25, 1-13
8+
!
9+
!
10+
program caroll_driver
11+
use stiff3_solver, only: wp => stiff3_wp, stiff3
12+
use caroll_test_problems
13+
14+
implicit none
15+
16+
integer :: n
17+
real(wp) :: x, xend, h0, eps, gerr
18+
real(wp) :: w(4), y(4)
19+
integer :: stats(6)
20+
21+
w = 1
22+
eps = 0.001_wp
23+
24+
!
25+
! P1
26+
!
27+
n = 3
28+
x = 0; xend = 100
29+
y(1:3) = [0,0,0]
30+
h0 = eps/20
31+
gerr = 0
32+
33+
call stiff3(n,p1_fun,x,y,xend,p1_jac,h0,eps,w,&
34+
stats=stats,solout=p1_out)
35+
call print_stats("P1",stats)
36+
print *, "Global error: ", gerr, gerr < eps
37+
38+
!
39+
! P2
40+
!
41+
42+
n = 3
43+
x = 0; xend = 20
44+
y(1:3) = [1,1,0]
45+
h0 = eps/20
46+
gerr = 0
47+
48+
call stiff3(n,p2_fun,x,y,xend,p2_jac,h0,eps,w,&
49+
stats=stats,solout=p2_out)
50+
call print_stats("P2",stats)
51+
print *, "Global error: ", gerr, gerr < eps
52+
53+
!
54+
! P3
55+
!
56+
57+
n = 3
58+
x = 0; xend = 25
59+
y(1:3) = [25498.0_wp/1500.0_wp,-16499.0_wp/1500.0_wp,x]
60+
h0 = eps/20
61+
gerr = 0
62+
63+
call stiff3(n,p3_fun,x,y,xend,p3_jac,h0,eps,w,&
64+
stats=stats,solout=p3_out)
65+
call print_stats("P3",stats)
66+
print *, "Global error: ", gerr, gerr < eps
67+
68+
!
69+
! P4
70+
!
71+
n = 3
72+
x = 0; xend = 20
73+
y(1:3) = [real(wp) :: 2,1,0]
74+
h0 = eps/20
75+
gerr = 0
76+
77+
call stiff3(n,p4_fun,x,y,xend,p4_jac,h0,eps,w,&
78+
stats=stats,solout=p4_out)
79+
call print_stats("P4",stats)
80+
print *, "Global error: ", gerr, gerr < eps
81+
82+
!
83+
! P5
84+
!
85+
n = 4
86+
x = 0; xend = 20
87+
y(1:4) = -1.0_wp
88+
h0 = eps/20
89+
gerr = 0
90+
91+
call stiff3(n,p5_fun,x,y,xend,p5_jac,h0,eps,w,&
92+
stats=stats,solout=p5_out)
93+
call print_stats("P5",stats)
94+
print *, "Global error: ", gerr, gerr < eps
95+
96+
!
97+
! P6 (Robertson)
98+
!
99+
n = 3
100+
x = 0; xend = 20
101+
y(1:3) = [1,0,0]
102+
h0 = eps/20
103+
call stiff3(n,p6_fun,x,y,xend,p6_jac,h0,eps,w,&
104+
stats=stats)
105+
call print_stats("P6",stats)
106+
107+
!
108+
! P7
109+
!
110+
n = 3
111+
x = 0; xend = 400
112+
y(1:3) = [0,0,0]
113+
h0 = eps/20
114+
call stiff3(n,p7_fun,x,y,xend,p7_jac,h0,eps,w,&
115+
stats=stats)
116+
call print_stats("P7",stats)
117+
118+
!
119+
! P8
120+
!
121+
n = 2
122+
x = 0; xend = 20
123+
y(1:2) = [1,0]
124+
h0 = eps/20
125+
call stiff3(n,p8_fun,x,y,xend,p8_jac,h0,eps,w,&
126+
stats=stats)
127+
call print_stats("P8",stats)
128+
129+
!
130+
! P9
131+
!
132+
n = 2
133+
x = 0; xend = 5
134+
y(1:2) = [10000.0_wp/10004.0_wp,1.0_wp]
135+
h0 = eps/20
136+
call stiff3(n,p9_fun,x,y,xend,p9_jac,h0,eps,w,&
137+
stats=stats)
138+
call print_stats("P9",stats)
139+
140+
!
141+
! P10
142+
!
143+
n = 2
144+
x = 0; xend = 20
145+
y(1:2) = [5,5]
146+
h0 = eps/20
147+
call stiff3(n,p10_fun,x,y,xend,p10_jac,h0,eps,w,&
148+
stats=stats)
149+
call print_stats("P10",stats)
150+
151+
152+
! ---------------------
153+
154+
!
155+
! P7, Convergence study with respect to step size
156+
! See Table 18, page 12, in Caroll (1989)
157+
!
158+
print '(A)', "P7 Convergence study"
159+
block
160+
integer :: m
161+
real(wp) :: x, xend, h0, y(3)
162+
n = 3
163+
eps = 0.1_wp
164+
do m = 3, 7
165+
x = 0; xend = 400
166+
y(1:3) = [0,0,0]
167+
h0 = 1.0_wp/real(2**m,wp)
168+
call stiff3(n,p7_fun,x,y,xend,p7_jac,h0,eps,w,&
169+
stats=stats,hmax=(h0))
170+
print *, y(1:3), stats
171+
end do
172+
end block
173+
174+
contains
175+
176+
subroutine print_stats(problem,stats)
177+
character(len=*), intent(in) :: problem
178+
integer, intent(in) :: stats(6)
179+
180+
print '(A)', problem
181+
print '(A,3(I0,1X))', 'accepted rejected nfev: ', stats(1), stats(2), stats(3)
182+
print '(A,3(I0,1X))', 'njev nlu nsol: ', stats(4), stats(5), stats(6)
183+
184+
end subroutine
185+
186+
subroutine p1_out(nr,told,t,y,ih,qa,irtrn)
187+
integer, intent(in) :: nr
188+
real(wp), intent(in) :: told
189+
real(wp), intent(in) :: t
190+
real(wp), intent(in) :: y(:)
191+
integer, intent(in) :: ih
192+
real(wp), intent(in) :: qa
193+
integer, intent(inout) :: irtrn
194+
195+
associate(sol => p1_sol(t))
196+
gerr = max(gerr,maxval(abs(y - sol),dim=1))
197+
end associate
198+
199+
end subroutine
200+
201+
subroutine p2_out(nr,told,t,y,ih,qa,irtrn)
202+
integer, intent(in) :: nr
203+
real(wp), intent(in) :: told
204+
real(wp), intent(in) :: t
205+
real(wp), intent(in) :: y(:)
206+
integer, intent(in) :: ih
207+
real(wp), intent(in) :: qa
208+
integer, intent(inout) :: irtrn
209+
210+
associate(sol => p2_sol(t))
211+
gerr = max(gerr,maxval(abs(y - sol),dim=1))
212+
end associate
213+
214+
end subroutine
215+
216+
subroutine p3_out(nr,told,t,y,ih,qa,irtrn)
217+
integer, intent(in) :: nr
218+
real(wp), intent(in) :: told
219+
real(wp), intent(in) :: t
220+
real(wp), intent(in) :: y(:)
221+
integer, intent(in) :: ih
222+
real(wp), intent(in) :: qa
223+
integer, intent(inout) :: irtrn
224+
225+
associate(sol => p3_sol(t))
226+
gerr = max(gerr,maxval(abs(y - sol),dim=1))
227+
end associate
228+
229+
end subroutine
230+
231+
subroutine p4_out(nr,told,t,y,ih,qa,irtrn)
232+
integer, intent(in) :: nr
233+
real(wp), intent(in) :: told
234+
real(wp), intent(in) :: t
235+
real(wp), intent(in) :: y(:)
236+
integer, intent(in) :: ih
237+
real(wp), intent(in) :: qa
238+
integer, intent(inout) :: irtrn
239+
240+
associate(sol => p4_sol(t))
241+
print *, t, abs(y - sol)
242+
gerr = max(gerr,maxval(abs(y - sol),dim=1))
243+
end associate
244+
245+
end subroutine
246+
247+
subroutine p5_out(nr,told,t,y,ih,qa,irtrn)
248+
integer, intent(in) :: nr
249+
real(wp), intent(in) :: told
250+
real(wp), intent(in) :: t
251+
real(wp), intent(in) :: y(:)
252+
integer, intent(in) :: ih
253+
real(wp), intent(in) :: qa
254+
integer, intent(inout) :: irtrn
255+
256+
associate(sol => p5_sol(t))
257+
gerr = max(gerr,maxval(abs(y - sol),dim=1))
258+
end associate
259+
260+
end subroutine
261+
end program

0 commit comments

Comments
 (0)