diff --git a/README.md b/README.md index 9fe1a6b..cfb856d 100644 --- a/README.md +++ b/README.md @@ -74,7 +74,7 @@ stiff3.git = "https://github.com/ivan-pi/stiff3" ### Running the examples -Six examples called `robertson`, `vanpol`, `lorenz`, `predator_prey`, `pendant_drop`, and `three_equation_system` are provided. They can be run with the command +Seven examples called `robertson`, `vanpol`, `lorenz`, `predator_prey`, `pendant_drop`, `three_equation_system`, and `oregonator` are provided. They can be run with the command ```sh fpm run --example @@ -82,6 +82,8 @@ fpm run --example With CMake, the compiled executables are placed in the `build/example/` directory and can be run directly or via CTest. +The `oregonator` example verifies a published reference solution at `t = 360` and writes `oregonator_work_precision.csv` together with an `oregonator_work_precision.svg` plot in the current working directory. + ## Usage Basic use of the solver is demonstrated using the [Van der Pol oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator): diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index dd2582a..c1a178d 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -1,4 +1,4 @@ -foreach(example robertson vanpol lorenz predator_prey pendant_drop three_equation_system) +foreach(example robertson vanpol lorenz predator_prey pendant_drop three_equation_system oregonator) add_executable(${example} ${example}.f90) target_link_libraries(${example} PRIVATE stiff3) add_test(NAME ${example} COMMAND ${example}) diff --git a/example/oregonator.f90 b/example/oregonator.f90 new file mode 100644 index 0000000..f2b74c9 --- /dev/null +++ b/example/oregonator.f90 @@ -0,0 +1,239 @@ +program oregonator + + use stiff3_solver, only: stiff3, wp => stiff3_wp + + implicit none + + integer, parameter :: n = 3 + integer, parameter :: ncases = 6 + real(wp), parameter :: y0(n) = [1.0_wp, 2.0_wp, 3.0_wp] + real(wp), parameter :: yref(n) = [ & + 1.000814870318523_wp, 1228.178521549917_wp, 132.0554942846706_wp ] + real(wp), parameter :: x0 = 0.0_wp + real(wp), parameter :: x1 = 360.0_wp + real(wp), parameter :: h0_initial = 1.0e-3_wp + real(wp), parameter :: verification_tol = 1.0e-8_wp + real(wp), parameter :: eps_values(ncases) = [ & + 1.0e-4_wp, 1.0e-5_wp, 1.0e-6_wp, 1.0e-7_wp, 1.0e-8_wp, 1.0e-9_wp ] + character(len=*), parameter :: csv_file = 'oregonator_work_precision.csv' + character(len=*), parameter :: svg_file = 'oregonator_work_precision.svg' + + integer :: i + integer :: verification_stats(6) + integer :: stats(6,ncases) + real(wp) :: y(n), w(n), errors(ncases), work(ncases), verification_error + + w = 1.0_wp + + call integrate_case(1.0e-8_wp, y, verification_stats) + verification_error = max_relative_error(y, yref) + if (verification_error > verification_tol) then + print '(A,ES12.4,A,ES12.4)', & + 'reference verification failed: max relative error=', verification_error, & + ' exceeds ', verification_tol + error stop 1 + end if + + print '(A,ES12.4)', 'verified reference solution with max relative error ', verification_error + print '(A,3(1X,ES24.16))', 'y(360) =', y + print '(A,6(I0,1X))', 'verification stats [nacc nrej nfev njev nlu nsol]: ', verification_stats + + do i = 1, ncases + call integrate_case(eps_values(i), y, stats(:,i)) + errors(i) = max_relative_error(y, yref) + work(i) = real(stats(3,i), wp) + end do + + call write_csv(csv_file, eps_values, errors, stats) + call write_svg(svg_file, eps_values, errors, work) + + print '(A)', 'work-precision data:' + print '(A)', ' eps error nfev njev nlu nsol' + do i = 1, ncases + print '(2(1X,ES10.2),4(1X,I6))', eps_values(i), errors(i), & + stats(3,i), stats(4,i), stats(5,i), stats(6,i) + end do + print '(A,1X,A)', 'wrote', csv_file + print '(A,1X,A)', 'wrote', svg_file + +contains + + subroutine integrate_case(eps, y, stats) + real(wp), intent(in) :: eps + real(wp), intent(out) :: y(n) + integer, intent(out) :: stats(6) + + real(wp) :: h0 + + y = y0 + h0 = h0_initial + call stiff3(n, fun, x0, y, x1, jac, h0, eps, w, stats=stats) + end subroutine + + function max_relative_error(y, y_reference) result(err) + real(wp), intent(in) :: y(n), y_reference(n) + real(wp) :: err + + err = maxval(abs(y - y_reference)/max(1.0_wp, abs(y_reference))) + end function + + subroutine write_csv(filename, eps, err, stats) + character(len=*), intent(in) :: filename + real(wp), intent(in) :: eps(:), err(:) + integer, intent(in) :: stats(:,:) + + integer :: i, ios, unit + + open(newunit=unit, file=filename, status='replace', action='write', iostat=ios) + if (ios /= 0) error stop 'failed to open CSV output file' + + write(unit,'(A)') 'eps,max_relative_error,nfev,njev,nlu,nsol,nacc,nrej' + do i = 1, size(eps) + write(unit,'(ES24.16,'','',ES24.16,'','',I0,'','',I0,'','',I0,'','',I0,'','',I0,'','',I0)') & + eps(i), err(i), stats(3,i), stats(4,i), stats(5,i), stats(6,i), stats(1,i), stats(2,i) + end do + + close(unit) + end subroutine + + subroutine write_svg(filename, eps, err, work) + character(len=*), intent(in) :: filename + real(wp), intent(in) :: eps(:), err(:), work(:) + + integer, parameter :: width = 800 + integer, parameter :: height = 600 + real(wp), parameter :: left = 100.0_wp + real(wp), parameter :: right = 30.0_wp + real(wp), parameter :: top = 50.0_wp + real(wp), parameter :: bottom = 80.0_wp + + character(len=32) :: label + integer :: i, ios, unit + real(wp) :: log_work(size(work)), log_err(size(err)) + real(wp) :: xmin, xmax, ymin, ymax, xpad, ypad + real(wp) :: x1p, x2p, y1p, y2p, xtick, ytick, tick_value + + log_work = log10(max(work, tiny(1.0_wp))) + log_err = log10(max(err, tiny(1.0_wp))) + + xmin = minval(log_work) + xmax = maxval(log_work) + ymin = minval(log_err) + ymax = maxval(log_err) + xpad = 0.05_wp*max(xmax - xmin, 1.0_wp) + ypad = 0.08_wp*max(ymax - ymin, 1.0_wp) + xmin = xmin - xpad + xmax = xmax + xpad + ymin = ymin - ypad + ymax = ymax + ypad + + open(newunit=unit, file=filename, status='replace', action='write', iostat=ios) + if (ios /= 0) error stop 'failed to open SVG output file' + + write(unit,'(A)') '' + write(unit,'(A,I0,A,I0,A)') '' + write(unit,'(A)') '' + write(unit,'(A)') & + 'Oregonator work-precision diagram' + write(unit,'(A,F0.3,A,F0.3,A,F0.3,A,F0.3,A)') & + '' + + do i = 0, 4 + xtick = left + real(i,wp)*(real(width,wp) - left - right)/4.0_wp + tick_value = 10.0_wp**(xmin + real(i,wp)*(xmax - xmin)/4.0_wp) + write(label,'(ES10.2)') tick_value + write(unit,'(A,F0.3,A,F0.3,A)') '' + write(unit,'(A,F0.3,A)') '' + write(unit,'(A)') trim(adjustl(label))//'' + end do + + do i = 0, 4 + ytick = real(height,wp) - bottom - real(i,wp)*(real(height,wp) - top - bottom)/4.0_wp + tick_value = 10.0_wp**(ymin + real(i,wp)*(ymax - ymin)/4.0_wp) + write(label,'(ES10.2)') tick_value + write(unit,'(A,F0.3,A,F0.3,A)') '' + write(unit,'(A,F0.3,A)') '' + write(unit,'(A)') trim(adjustl(label))//'' + end do + + write(unit,'(A)') & + 'rhs evaluations (nfev)' + write(unit,'(A)') & + '' // & + 'max relative error' + + do i = 1, size(work) - 1 + x1p = xcoord(log_work(i), xmin, xmax, width, left, right) + y1p = ycoord(log_err(i), ymin, ymax, height, top, bottom) + x2p = xcoord(log_work(i + 1), xmin, xmax, width, left, right) + y2p = ycoord(log_err(i + 1), ymin, ymax, height, top, bottom) + write(unit,'(A,F0.3,A,F0.3,A,F0.3,A,F0.3,A)') '' + end do + + do i = 1, size(work) + x1p = xcoord(log_work(i), xmin, xmax, width, left, right) + y1p = ycoord(log_err(i), ymin, ymax, height, top, bottom) + write(label,'(ES9.1)') eps(i) + write(unit,'(A,F0.3,A,F0.3,A)') '' + write(unit,'(A,F0.3,A,F0.3,A)') '' + write(unit,'(A)') 'eps='//trim(adjustl(label))//'' + end do + + write(unit,'(A)') '' + + close(unit) + + end subroutine + + function xcoord(value, xmin, xmax, width, left, right) result(xplot) + real(wp), intent(in) :: value, xmin, xmax, left, right + integer, intent(in) :: width + real(wp) :: xplot + + xplot = left + (value - xmin)*(real(width,wp) - left - right)/(xmax - xmin) + end function + + function ycoord(value, ymin, ymax, height, top, bottom) result(yplot) + real(wp), intent(in) :: value, ymin, ymax, top, bottom + integer, intent(in) :: height + real(wp) :: yplot + + yplot = real(height,wp) - bottom - & + (value - ymin)*(real(height,wp) - top - bottom)/(ymax - ymin) + end function + + subroutine fun(n, y, f) + integer, intent(in) :: n + real(wp), intent(in) :: y(n) + real(wp), intent(inout) :: f(n) + + f(1) = 77.27_wp*(y(2) + y(1)*(1.0_wp - 8.375e-6_wp*y(1) - y(2))) + f(2) = (y(3) - (1.0_wp + y(1))*y(2))/77.27_wp + f(3) = 0.161_wp*(y(1) - y(3)) + end subroutine + + subroutine jac(n, y, df) + integer, intent(in) :: n + real(wp), intent(in) :: y(n) + real(wp), intent(inout) :: df(n,n) + + df = 0.0_wp + df(1,1) = 77.27_wp*(1.0_wp - 2.0_wp*8.375e-6_wp*y(1) - y(2)) + df(1,2) = 77.27_wp*(1.0_wp - y(1)) + df(2,1) = -y(2)/77.27_wp + df(2,2) = -(1.0_wp + y(1))/77.27_wp + df(2,3) = 1.0_wp/77.27_wp + df(3,1) = 0.161_wp + df(3,3) = -0.161_wp + end subroutine + +end program