Skip to content

Commit 9c23db8

Browse files
Copilotivan-pi
Copilot
and
ivan-pi
authored
Add HIRES example with eps sweep output data (#46)
Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/19e7cc0c-039e-4101-9be5-cb259db8f357 Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>
1 parent 6fb3a96 commit 9c23db8

3 files changed

Lines changed: 159 additions & 2 deletions

File tree

‎README.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ stiff3.git = "https://github.com/ivan-pi/stiff3"
7474

7575
### Running the examples
7676

77-
Eight examples called `robertson`, `vanpol`, `lorenz`, `predator_prey`, `pendant_drop`, `three_equation_system`, `oregonator`, and `e5_chemical` are provided. They can be run with the command
77+
Nine examples called `robertson`, `vanpol`, `lorenz`, `predator_prey`, `pendant_drop`, `three_equation_system`, `oregonator`, `e5_chemical`, and `hires` are provided. They can be run with the command
7878

7979
```sh
8080
fpm run --example <name>
@@ -86,6 +86,8 @@ The `oregonator` example verifies a published reference solution at `t = 360` an
8686

8787
The `e5_chemical` example verifies the Datta/Bari E5 reference solution at `t = 1.0e13` and writes `e5_chemical_work_precision.csv` together with an `e5_chemical_work_precision.svg` plot of CPU time vs error tolerance `eps` (logarithmic y-axis).
8888

89+
The `hires` example verifies the Bari HIRES reference solution at `t = 321.8122` and writes `hires_work_precision.dat` (columns: `eps cpu_time_seconds max_relative_error nfev njev nlu nsol nacc nrej`) for plotting CPU time vs error tolerance `eps` with gnuplot.
90+
8991
## Usage
9092

9193
Basic use of the solver is demonstrated using the [Van der Pol oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator):

‎example/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
foreach(example robertson vanpol lorenz predator_prey pendant_drop three_equation_system oregonator e5_chemical)
1+
foreach(example robertson vanpol lorenz predator_prey pendant_drop three_equation_system oregonator e5_chemical hires)
22
add_executable(${example} ${example}.f90)
33
target_link_libraries(${example} PRIVATE stiff3)
44
add_test(NAME ${example} COMMAND ${example})

‎example/hires.f90‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
program hires
2+
3+
use stiff3_solver, only: stiff3, wp => stiff3_wp
4+
5+
implicit none
6+
7+
integer, parameter :: n = 8
8+
integer, parameter :: ncases = 6
9+
real(wp), parameter :: x0 = 0.0_wp
10+
real(wp), parameter :: x1 = 321.8122_wp
11+
real(wp), parameter :: h0_initial = 1.0e-4_wp
12+
real(wp), parameter :: verification_tol = 1.0e-8_wp
13+
real(wp), parameter :: eps_values(ncases) = [ &
14+
1.0e-4_wp, 1.0e-5_wp, 1.0e-6_wp, 1.0e-7_wp, 1.0e-8_wp, 1.0e-9_wp ]
15+
real(wp), parameter :: y0(n) = [ &
16+
1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 5.7e-3_wp ]
17+
real(wp), parameter :: yref(n) = [ &
18+
0.7371312573325668e-3_wp, &
19+
0.1442485726316185e-3_wp, &
20+
0.5888729740967575e-4_wp, &
21+
0.1175651343283149e-2_wp, &
22+
0.2386356198831331e-2_wp, &
23+
0.6238968252742796e-2_wp, &
24+
0.2849998395185769e-2_wp, &
25+
0.2850001604814231e-2_wp ]
26+
character(len=*), parameter :: output_file = 'hires_work_precision.dat'
27+
28+
integer :: i
29+
integer :: verification_stats(6)
30+
integer :: stats(6,ncases)
31+
real(wp) :: y(n), w(n), errors(ncases), cpu_times(ncases), verification_error
32+
33+
w = 1.0_wp
34+
35+
call integrate_case(1.0e-9_wp, y, verification_stats, verification_error)
36+
if (verification_error > verification_tol) then
37+
print '(A,ES12.4,A,ES12.4)', &
38+
'reference verification failed: max relative error=', verification_error, &
39+
' exceeds ', verification_tol
40+
error stop 1
41+
end if
42+
43+
print '(A,ES12.4)', 'verified reference solution with max relative error ', verification_error
44+
print '(A,8(1X,ES24.16))', 'y(321.8122) =', y
45+
print '(A,6(I0,1X))', 'verification stats [nacc nrej nfev njev nlu nsol]: ', verification_stats
46+
47+
do i = 1, ncases
48+
call integrate_case(eps_values(i), y, stats(:,i), errors(i), cpu_times(i))
49+
end do
50+
51+
call write_gnuplot_data(output_file, eps_values, cpu_times, errors, stats)
52+
53+
print '(A)', 'work-precision data (CPU time vs eps):'
54+
print '(A)', ' eps cpu[s] error nfev njev nlu nsol'
55+
do i = 1, ncases
56+
print '(1X,ES10.2,1X,ES12.4,1X,ES10.2,4(1X,I6))', eps_values(i), cpu_times(i), errors(i), &
57+
stats(3,i), stats(4,i), stats(5,i), stats(6,i)
58+
end do
59+
print '(A,1X,A)', 'wrote', output_file
60+
61+
contains
62+
63+
subroutine integrate_case(eps, y, stats, err, cpu_time_seconds)
64+
real(wp), intent(in) :: eps
65+
real(wp), intent(out) :: y(n)
66+
integer, intent(out) :: stats(6)
67+
real(wp), intent(out), optional :: err, cpu_time_seconds
68+
69+
real(wp), parameter :: min_cpu_time = 1.0e-12_wp
70+
real(wp) :: h0, cpu_start, cpu_end
71+
72+
y = y0
73+
h0 = h0_initial
74+
call cpu_time(cpu_start)
75+
call stiff3(n, fun, x0, y, x1, jac, h0, eps, w, stats=stats)
76+
call cpu_time(cpu_end)
77+
if (present(err)) err = max_relative_error(y, yref)
78+
if (present(cpu_time_seconds)) cpu_time_seconds = max(min_cpu_time, cpu_end - cpu_start)
79+
end subroutine
80+
81+
function max_relative_error(y, y_reference) result(err)
82+
real(wp), intent(in) :: y(n), y_reference(n)
83+
real(wp) :: err
84+
85+
err = maxval(abs(y - y_reference)/max(1.0_wp, abs(y_reference)))
86+
end function
87+
88+
subroutine write_gnuplot_data(filename, eps, cpu, err, stats)
89+
character(len=*), intent(in) :: filename
90+
real(wp), intent(in) :: eps(:), cpu(:), err(:)
91+
integer, intent(in) :: stats(:,:)
92+
93+
integer :: i, ios, unit
94+
95+
open(newunit=unit, file=filename, status='replace', action='write', iostat=ios)
96+
if (ios /= 0) error stop 'failed to open HIRES output file'
97+
98+
write(unit,'(A)') '# eps cpu_time_seconds max_relative_error nfev njev nlu nsol nacc nrej'
99+
do i = 1, size(eps)
100+
write(unit,'(ES24.16,1X,ES24.16,1X,ES24.16,1X,I0,1X,I0,1X,I0,1X,I0,1X,I0,1X,I0)') &
101+
eps(i), cpu(i), err(i), stats(3,i), stats(4,i), stats(5,i), stats(6,i), stats(1,i), stats(2,i)
102+
end do
103+
104+
close(unit)
105+
end subroutine
106+
107+
subroutine fun(n, y, f)
108+
integer, intent(in) :: n
109+
real(wp), intent(in) :: y(n)
110+
real(wp), intent(inout) :: f(n)
111+
112+
f(1) = -1.71_wp*y(1) + 0.43_wp*y(2) + 8.32_wp*y(3) + 0.0007_wp
113+
f(2) = 1.71_wp*y(1) - 8.75_wp*y(2)
114+
f(3) = -10.03_wp*y(3) + 0.43_wp*y(4) + 0.035_wp*y(5)
115+
f(4) = 8.32_wp*y(2) + 1.71_wp*y(3) - 1.12_wp*y(4)
116+
f(5) = -1.745_wp*y(5) + 0.43_wp*(y(6) + y(7))
117+
f(6) = -280.0_wp*y(6)*y(8) + 0.69_wp*y(4) + 1.71_wp*y(5) - 0.43_wp*y(6) + 0.69_wp*y(7)
118+
f(7) = 280.0_wp*y(6)*y(8) - 1.81_wp*y(7)
119+
f(8) = -f(7)
120+
end subroutine
121+
122+
subroutine jac(n, y, df)
123+
integer, intent(in) :: n
124+
real(wp), intent(in) :: y(n)
125+
real(wp), intent(inout) :: df(n,n)
126+
127+
df = 0.0_wp
128+
df(1,1) = -1.71_wp
129+
df(1,2) = 0.43_wp
130+
df(1,3) = 8.32_wp
131+
df(2,1) = 1.71_wp
132+
df(2,2) = -8.75_wp
133+
df(3,3) = -10.03_wp
134+
df(3,4) = 0.43_wp
135+
df(3,5) = 0.035_wp
136+
df(4,2) = 8.32_wp
137+
df(4,3) = 1.71_wp
138+
df(4,4) = -1.12_wp
139+
df(5,5) = -1.745_wp
140+
df(5,6) = 0.43_wp
141+
df(5,7) = 0.43_wp
142+
df(6,4) = 0.69_wp
143+
df(6,5) = 1.71_wp
144+
df(6,6) = -280.0_wp*y(8) - 0.43_wp
145+
df(6,7) = 0.69_wp
146+
df(6,8) = -280.0_wp*y(6)
147+
df(7,6) = 280.0_wp*y(8)
148+
df(7,7) = -1.81_wp
149+
df(7,8) = 280.0_wp*y(6)
150+
df(8,6) = -280.0_wp*y(8)
151+
df(8,7) = 1.81_wp
152+
df(8,8) = -280.0_wp*y(6)
153+
end subroutine
154+
155+
end program

0 commit comments

Comments
 (0)