From 217a1fb2b133e54d6d4364ec9adee6db61e9964b Mon Sep 17 00:00:00 2001 From: Niels Claes Date: Mon, 17 Apr 2023 14:31:30 +0200 Subject: [PATCH] make RK45 integration conform with function behaviour --- src/mod_integration.f08 | 405 ++++++----------------- src/physics/mod_solar_atmosphere.f08 | 38 ++- tests/unit_tests/mod_test_integration.pf | 249 +++++--------- 3 files changed, 219 insertions(+), 473 deletions(-) diff --git a/src/mod_integration.f08 b/src/mod_integration.f08 index dda492a8..a38850df 100644 --- a/src/mod_integration.f08 +++ b/src/mod_integration.f08 @@ -6,338 +6,149 @@ !! $$ y'(x) = A(x)y(x) + B(x) $$ !! These are solved using a fifth-order Runge-Kutta method. module mod_integration - use mod_global_variables, only: dp, dp_LIMIT - use mod_logging, only: logger, str + use mod_global_variables, only: dp + use mod_logging, only: logger, str, exp_fmt implicit none private - public :: integrate_ode_rk + interface + real(dp) function func(x) + use mod_global_variables, only: dp + real(dp), intent(in) :: x + end function func + end interface -contains + public :: integrate_ode_rk45 +contains - !> Integrates a first order differential equation of the form - !! $$ y'(x) = A(x)y(x) + B(x) $$ using a fifth-order Runge-Kutta - !! method. The argument nbpoints determines the stepsize through - !! $$ dh = \frac{xvalues(N) - xvalues(1)}{nbpoints} $$ - !! If the arrays \(A(x), B(x)\) are not of size nbpoints, then these are - !! interpolated to that resolution. The differential equation is then integrated, - !! the solution will also be of size nbpoints and can be downsampled using - !! the appropriate subroutine. - !! If desired the optional argument dyvalues can be provided, - !! which contains the (numerical) derivative of y. - subroutine integrate_ode_rk( & - xvalues, & - axvalues, & - bxvalues, & - nbpoints, & - yinit, & - yvalues, & - adaptive, & - epsilon, & - max_step_change, & - min_dh_size, & - max_dh_size, & - max_iter_per_step, & - new_xvalues & + !> Numerically integrates the differential equation + !! $$ y'(x) = A(x)y(x) + B(x) $$ using a fifth-order Runge-Kutta method. + !! The functions A(x) and B(x) are passed as arguments and should be conform to the + !! given interface, that is, these should be `real(dp)` functions which take a single + !! `real(dp), intent(in)` argument. + !! The integration is performed on the interval [x0, x1] with a stepsize of + !! `dh = (x1 - x0) / (nbpoints - 1)`. + subroutine integrate_ode_rk45( & + x0, x1, ax_func, bx_func, nbpoints, yinit, yvalues, xvalues & ) - use mod_interpolation, only: interpolate_table - - !> array of x-values - real(dp), intent(in) :: xvalues(:) - !> term \(A(x)\) - real(dp), intent(in) :: axvalues(:) - !> term \(B(x)\) - real(dp), intent(in) :: bxvalues(:) - !> desired resolution of coefficient arrays axvalues and bxvalues - integer, intent(in) :: nbpoints - !> initial value for \(y\) at left edge (xvalues(1)) - real(dp), intent(in) :: yinit - !> array of y-values, size will be nbpoints if adaptive=.false., - !! size will vary if adaptive=.true. + !> start of x-interval + real(dp), intent(in) :: x0 + !> end of x-interval + real(dp), intent(in) :: x1 + !> function to calculate A(x) + procedure(func) :: ax_func + !> function to calculate B(x) + procedure(func) :: bx_func + !> number of points, determines stepsize + integer, intent(in) :: nbpoints + !> initial value of y + real(dp), intent(in) :: yinit + !> integrated values of y real(dp), intent(out), allocatable :: yvalues(:) - !> if .true. an adaptive stepsize will be used - logical, intent(in), optional :: adaptive - !> maximum truncation error, determines adaptive step - real(dp), intent(in), optional :: epsilon - !> maximum factor by which dh can change if step is adaptive - real(dp), intent(in), optional :: max_step_change - !> limit on the minimal stepsize - real(dp), intent(in), optional :: min_dh_size - !> limit on the maximal stepsize - real(dp), intent(in), optional :: max_dh_size - !> limit on maximum number of iterations for a single step - integer, intent(in), optional :: max_iter_per_step - !> array of x-values corresponding to yvalues - real(dp), intent(out), allocatable, optional :: new_xvalues(:) - - integer :: i, re_iter, maxiter - real(dp) :: ax_interp(nbpoints), bx_interp(nbpoints), x_interp(nbpoints) - real(dp) :: xi, xend, yi, dh, dh_new, dh_min, dh_max - real(dp) :: ysolrk4, ysolrk5, err, tol, max_dh_fact - logical :: use_adaptive_stepping - real(dp), allocatable :: xtemp(:), ytemp(:) - - ! check if A(x) needs to be resampled - if (needs_resampling(size(axvalues), nbpoints, "A(x)")) then - call interpolate_table(nbpoints, xvalues, axvalues, x_interp, ax_interp) - else - ax_interp = axvalues - x_interp = xvalues - end if - ! check if B(x) needs to be resampled - if (needs_resampling(size(bxvalues), nbpoints, "B(x)")) then - call interpolate_table(nbpoints, xvalues, bxvalues, x_interp, bx_interp) - else - bx_interp = bxvalues - x_interp = xvalues - end if - - ! determine starting stepsize - dh = x_interp(2) - x_interp(1) - ! set max dh change - max_dh_fact = 1.5d0 - if (present(max_step_change)) then - max_dh_fact = max_step_change - end if - ! check adaptive stepping tolerance - tol = 1.0d-4 - if (present(epsilon)) then - tol = epsilon - end if - ! check min/max stepsizes - dh_min = 1.0d-10 - dh_max = dh - maxiter = 50 - if (present(min_dh_size)) then - dh_min = min_dh_size - end if - if (present(max_dh_size)) then - dh_max = max_dh_size - end if - if (present(max_iter_per_step)) then - maxiter = max_iter_per_step - end if - ! set adaptive stepping - use_adaptive_stepping = .false. - if (present(adaptive)) then - use_adaptive_stepping = adaptive - if (use_adaptive_stepping) then - call logger%info( & - "integrating using adaptive stepping, max dh change = "// str(max_dh_fact) & - ) - else - call logger%info("integrating using regular stepping, dh = " // str(dh)) - end if - end if - - ! allocate temporary arrays, this should be more than sufficient - if (use_adaptive_stepping) then - allocate(xtemp(10 * nbpoints), ytemp(10 * nbpoints)) - else - allocate(xtemp(nbpoints), ytemp(nbpoints)) - end if - - ! do integration - i = 1 - xi = x_interp(1) - xend = x_interp(nbpoints) - xtemp(1) = xi - ytemp(1) = yinit - re_iter = 0 - do while (xi < xend) - dh = min(dh, xend - xi) - yi = ytemp(i) - - ! sanity check, for non-adaptive stepping we sometimes we have xend - xi - ! equal to about 1e-13 at the last iteration i == nbpoints, and the above part - ! selects that as dh near the end. This results in an "additional" iteration - ! with a bogus dh due to numerical round-off, and a final array of length - ! nbpoints + 1. Hence, if we're not using adaptive stepping and we encounter - ! something like this, modify dh so that it nicely reaches xend next iteration. - if ( & - .not. use_adaptive_stepping & - .and. i == nbpoints - 1 & - ) then - dh = xend - xi - end if - - call rk45(xi, dh, x_interp, ax_interp, bx_interp, yi, ysolrk4, ysolrk5) - err = abs(ysolrk5 - ysolrk4) / dh - - if (use_adaptive_stepping) then - ! calculate truncation error and new stepsize - dh_new = 0.9d0 * dh * (tol / err) ** 0.2d0 - ! make sure stepsize does not change too much - if (dh_new / dh > max_dh_fact) then - dh = dh * max_dh_fact ! stepsize was increased too much - else if (dh / dh_new > max_dh_fact) then - dh = dh / max_dh_fact ! stepsize was decreased too much - else - dh = dh_new ! stepsize is fine, replace - end if - ! make sure stepsize is within bounds - if (dh < dh_min) then - dh = dh_min - else if (dh > dh_max) then - dh = dh_max - end if - ! if error is not within tolerance, retake this step with new dh - ! continue if iterations is higher than allowed per step - if (err > tol .and. re_iter < maxiter) then - re_iter = re_iter + 1 - cycle - end if - re_iter = 0 - ! also handle tiny dh for adaptive stepping - if ((xend - xi - dh) < dh_min) then - dh = xend - xi - end if - end if - - ! we are either within tolerance or without adaptive stepping, take next step - i = i + 1 + !> x-values corresponding to integrated y-values + real(dp), intent(out), allocatable, optional :: xvalues(:) + + real(dp) :: xi, yi, dh + real(dp) :: ysolrk4, ysolrk5 + real(dp), allocatable :: xivalues(:) + integer :: i + + dh = (x1 - x0) / (nbpoints - 1) + call logger%info( & + "integrating from " // str(x0) // " to " // str(x1) & + // " with dh = " // str(dh, fmt=exp_fmt) & + ) + xi = x0 + yi = yinit + allocate(yvalues(nbpoints), xivalues(nbpoints)) + xivalues(1) = xi + yvalues(1) = yi + do i = 2, nbpoints + call rk45(xi, dh, ax_func, bx_func, yi, ysolrk4, ysolrk5) + yi = ysolrk5 xi = xi + dh - xtemp(i) = xi - ytemp(i) = ysolrk5 - ! LCOV_EXCL_START - if (mod(i, 10000) == 0) then - call logger%info( & - "steps: " // str(i) // " | xi: " // str(xi) // " | dh: " & - // str(dh, fmt="e20.5") & - ) - end if - ! LCOV_EXCL_STOP + yvalues(i) = yi + xivalues(i) = xi end do - ! allocate final arrays - allocate(yvalues(i)) - yvalues = ytemp(:i) - deallocate(ytemp) - if (present(new_xvalues)) then - allocate(new_xvalues(i)) - new_xvalues = xtemp(:i) - deallocate(xtemp) - end if - end subroutine integrate_ode_rk + if (present(xvalues)) xvalues = xivalues + deallocate(xivalues) + end subroutine integrate_ode_rk45 !> Calculates the Runge-Kutta coefficients and calculates the fourth and fifth !! order solutions for step i+1 based on the values at step i. - subroutine rk45(xi, dh, xarray, axarray, bxarray, yi, ysolrk4, ysolrk5) - use mod_interpolation, only: lookup_table_value - !> current x value in iteration - real(dp), intent(in) :: xi + subroutine rk45(xi, dh, ax_func, bx_func, yi, ysolrk4, ysolrk5) + !> current x value + real(dp), intent(in) :: xi !> current step size - real(dp), intent(in) :: dh - !> array of x values - real(dp), intent(in) :: xarray(:) - !> array of A(x) values - real(dp), intent(in) :: axarray(:) - !> array of B(x) values - real(dp), intent(in) :: bxarray(:) - !> current y value in iteration - real(dp), intent(in) :: yi - !> fourth-order solution + real(dp), intent(in) :: dh + !> function to calculate A(x) + procedure(func) :: ax_func + !> function to calculate B(x) + procedure(func) :: bx_func + !> current y value + real(dp), intent(in) :: yi + !> fourth order solution real(dp), intent(out) :: ysolrk4 - !> fifth-order solution + !> fifth order solution real(dp), intent(out) :: ysolrk5 - real(dp) :: axi, bxi, yvalrk - real(dp) :: rkf1, rkf2, rkf3, rkf4, rkf5, rkf6 + real(dp) :: rkf1, rkf2, rkf3, rkf4, rkf5, rkf6 + real(dp) :: xvalrk, yvalrk ! first step - axi = lookup_table_value(xi, xarray, axarray) - bxi = lookup_table_value(xi, xarray, bxarray) - yvalrk = yi - rkf1 = dh * (axi * yvalrk + bxi) + rkf1 = dh * (ax_func(xi) * yi + bx_func(xi)) ! second step - axi = lookup_table_value(xi + 0.25d0 * dh, xarray, axarray) - bxi = lookup_table_value(xi + 0.25d0 * dh, xarray, bxarray) - yvalrk = yi + 0.25d0 * rkf1 - rkf2 = dh * (axi * yvalrk + bxi) + xvalrk = xi + 0.25_dp * dh + yvalrk = yi + 0.25_dp * rkf1 + rkf2 = dh * (ax_func(xvalrk) * yvalrk + bx_func(xvalrk)) ! third step - axi = lookup_table_value(xi + 3.0d0 * dh / 8.0d0, xarray, axarray) - bxi = lookup_table_value(xi + 3.0d0 * dh / 8.0d0, xarray, bxarray) - yvalrk = yi + (3.0d0 * rkf1 + 9.0d0 * rkf2) / 32.0d0 - rkf3 = dh * (axi * yvalrk + bxi) + xvalrk = xi + 3.0_dp * dh / 8.0_dp + yvalrk = yi + (3.0_dp * rkf1 + 9.0_dp * rkf2) / 32.0_dp + rkf3 = dh * (ax_func(xvalrk) * yvalrk + bx_func(xvalrk)) ! fourth step - axi = lookup_table_value(xi + 12.0d0 * dh / 13.0d0, xarray, axarray) - bxi = lookup_table_value(xi + 12.0d0 * dh / 13.0d0, xarray, bxarray) - yvalrk = yi + (1932.0d0 * rkf1 - 7200.0d0 * rkf2 + 7296.0d0 * rkf3) / 2197.0d0 - rkf4 = dh * (axi * yvalrk + bxi) + xvalrk = xi + 12.0_dp * dh / 13.0_dp + yvalrk = yi + (1932.0_dp * rkf1 - 7200.0_dp * rkf2 + 7296.0_dp * rkf3) / 2197.0_dp + rkf4 = dh * (ax_func(xvalrk) * yvalrk + bx_func(xvalrk)) ! fifth step - axi = lookup_table_value(xi + dh, xarray, axarray) - bxi = lookup_table_value(xi + dh, xarray, bxarray) - yvalrk = ( & - yi & - + (439.0d0 / 216.0d0) * rkf1 & - - 8.0d0 * rkf2 & - + (3680.0d0 / 513.0d0) * rkf3 & - - (845.0d0 / 4104.0d0) * rkf4 & + xvalrk = xi + dh + yvalrk = yi + ( & + 439.0_dp * rkf1 / 216.0_dp & + - 8.0_dp * rkf2 & + + 3680.0_dp * rkf3 / 513.0_dp & + - 845.0_dp * rkf4 / 4104.0_dp & ) - rkf5 = dh * (axi * yvalrk + bxi) + rkf5 = dh * (ax_func(xvalrk) * yvalrk + bx_func(xvalrk)) ! sixth step - axi = lookup_table_value(xi + 0.5d0 * dh, xarray, axarray) - bxi = lookup_table_value(xi + 0.5d0 * dh, xarray, bxarray) - yvalrk = ( & - yi & - - (8.0d0 / 27.0d0) * rkf1 & - + 2.0d0 * rkf2 & - - (3544.0d0 / 2565.0d0) * rkf3 & - + (1859.0d0 / 4104.0d0) * rkf4 & - - (11.0d0 / 40.0d0) * rkf5 & + xvalrk = xi + 0.5_dp * dh + yvalrk = yi + ( & + - 8.0_dp * rkf1 / 27.0_dp & + + 2.0_dp * rkf2 & + - 3544.0_dp * rkf3 / 2565.0_dp & + + 1859.0_dp * rkf4 / 4104.0_dp & + - 11.0_dp * rkf5 / 40.0_dp & ) - rkf6 = dh * (axi * yvalrk + bxi) + rkf6 = dh * (ax_func(xvalrk) * yvalrk + bx_func(xvalrk)) ! fourth order solution - ysolrk4 = ( & - yi & - + (25.0d0 / 216.0d0) * rkf1 & - + (1408.0d0 / 2565.0d0) * rkf3 & - + (2197.0d0 / 4104.0d0) * rkf4 & - - (1.0d0 / 5.0d0) * rkf5 & + ysolrk4 = yi + ( & + 25.0_dp * rkf1 / 216.0_dp & + + 1408.0_dp * rkf3 / 2565.0_dp & + + 2197.0_dp * rkf4 / 4104.0_dp & + - rkf5 / 5.0_dp & ) ! fifth order solution - ysolrk5 = ( & - yi & - + (16.0d0 / 135.0d0) * rkf1 & - + (6656.0d0 / 12825.0d0) * rkf3 & - + (28561.0d0 / 56430.0d0) * rkf4 & - - (9.0d0 / 50.0d0) * rkf5 & - + (2.0d0 / 55.0d0) * rkf6 & + ysolrk5 = yi + ( & + 16.0_dp * rkf1 / 135.0_dp & + + 6656.0_dp * rkf3 / 12825.0_dp & + + 28561.0_dp * rkf4 / 56430.0_dp & + - 9.0_dp * rkf5 / 50.0_dp & + + 2.0_dp * rkf6 / 55.0_dp & ) end subroutine rk45 - !> Checks if an array needs resampling. - function needs_resampling(base, target, array_name) result(resample) - !> size of base array - integer, intent(in) :: base - !> size of target array - integer, intent(in) :: target - !> name of array (printed to console) - character(len=*), intent(in) :: array_name - !> is .true. if arrays differ in size - logical :: resample - - ! LCOV_EXCL_START - if (target < base) then - call logger%warning( & - "resampling: target #points is less than size of input arrays!" & - ) - end if - ! LCOV_EXCL_STOP - - if (base == target) then - resample = .false. - else - call logger%info( & - "ode integrator: resampling " // trim(adjustl(array_name)) // " (" & - // str(base) // " -> " // str(target) // " points)" & - ) - resample = .true. - end if - end function needs_resampling end module mod_integration diff --git a/src/physics/mod_solar_atmosphere.f08 b/src/physics/mod_solar_atmosphere.f08 index 506a8ee3..21b6850b 100644 --- a/src/physics/mod_solar_atmosphere.f08 +++ b/src/physics/mod_solar_atmosphere.f08 @@ -52,8 +52,7 @@ module mod_solar_atmosphere !! can all use the same result. !! @warning Throws an error if the geometry is not Cartesian. @endwarning subroutine set_solar_atmosphere(settings, background, physics, n_interp) - use mod_interpolation, only: lookup_table_value, get_numerical_derivative - use mod_integration, only: integrate_ode_rk + use mod_integration, only: integrate_ode_rk45 type(settings_t), intent(in) :: settings type(background_t), intent(inout) :: background @@ -61,9 +60,7 @@ subroutine set_solar_atmosphere(settings, background, physics, n_interp) !> points used for interpolation, defaults to 4000 if not present integer, intent(in), optional :: n_interp - integer :: i - real(dp) :: x, rhoinit - real(dp), allocatable :: axvalues(:), bxvalues(:) + real(dp) :: rhoinit unit_length = settings%units%get_unit_length() unit_time = settings%units%get_unit_time() @@ -85,25 +82,20 @@ subroutine set_solar_atmosphere(settings, background, physics, n_interp) ! curves are normalised on return call create_atmosphere_curves(settings) - ! fill ODE functions, use high resolution arrays - allocate(axvalues(nbpoints), bxvalues(nbpoints)) - do i = 1, nbpoints - x = h_interp(i) - axvalues(i) = -(dT_interp(i) + g0(x)) / T_interp(i) - bxvalues = -(B02() * dB02() + B03() * dB03()) / T_interp(i) - end do ! set initial density value (numberdensity = density in normalised units) rhoinit = nh_interp(1) - ! solve differential equation call logger%info( & "solving equilibrium ODE for density (" // str(nbpoints) // " points)" & ) - ! do integration - call integrate_ode_rk( & - h_interp, axvalues, bxvalues, nbpoints, rhoinit, rho_values, adaptive=.false. & + call integrate_ode_rk45( & + x0=h_interp(1), & + x1=h_interp(nbpoints), & + ax_func=ax_values, & + bx_func=bx_values, & + nbpoints=nbpoints, & + yinit=rhoinit, & + yvalues=rho_values & ) - ! these are no longer needed - deallocate(axvalues, bxvalues) call background%set_density_funcs(rho0_func=rho0, drho0_func=drho0) call background%set_temperature_funcs(T0_func=T0, dT0_func=dT0) @@ -162,6 +154,16 @@ subroutine create_atmosphere_curves(settings) end subroutine create_atmosphere_curves + real(dp) function ax_values(x) + real(dp), intent(in) :: x + ax_values = -(dT0(x) + g0(x)) / T0(x) + end function ax_values + + real(dp) function bx_values(x) + real(dp), intent(in) :: x + bx_values = -(B02() * dB02() + B03() * dB03()) / T0(x) + end function bx_values + real(dp) function B02() B02 = 0.0_dp end function B02 diff --git a/tests/unit_tests/mod_test_integration.pf b/tests/unit_tests/mod_test_integration.pf index 7558f86b..860f85da 100644 --- a/tests/unit_tests/mod_test_integration.pf +++ b/tests/unit_tests/mod_test_integration.pf @@ -1,200 +1,133 @@ module mod_test_integration use funit use mod_suite_utils - use mod_integration, only: integrate_ode_rk + use mod_integration, only: integrate_ode_rk45 implicit none - real(dp) :: yinit - real(dp), allocatable :: xvalues(:), yactual(:) - real(dp), allocatable :: ax(:), bx(:) - real(dp), allocatable :: xrkvalues(:), yexpect(:) + integer, parameter :: pts = 5000 + real(dp) :: yinit + real(dp), allocatable :: yactual(:), yexpect(:) + real(dp), allocatable :: xactual(:), xvalues(:) contains - subroutine reset_arrays(pts) - integer, intent(in) :: pts - - if (allocated(xvalues)) then - deallocate(xvalues) - end if - if (allocated(yactual)) then - deallocate(yactual) - end if - if (allocated(ax)) then - deallocate(ax) - end if - if (allocated(bx)) then - deallocate(bx) - end if - if (allocated(xrkvalues)) then - deallocate(xrkvalues) - end if - if (allocated(yexpect)) then - deallocate(yexpect) - end if - allocate(xvalues(pts)) + @before + subroutine init_test() + call reset_globals() allocate(yactual(pts)) - allocate(ax(pts)) - allocate(bx(pts)) - allocate(xrkvalues(pts)) allocate(yexpect(pts)) - end subroutine reset_arrays + allocate(xactual(pts)) + allocate(xvalues(pts)) + end subroutine init_test + + + @after + subroutine teardown_test() + if (allocated(yactual)) deallocate(yactual) + if (allocated(yexpect)) deallocate(yexpect) + if (allocated(xactual)) deallocate(xactual) + if (allocated(xvalues)) deallocate(xvalues) + end subroutine teardown_test @test subroutine test_integration_simple() call set_name("integration (dy=x, x=[1, 3], y0=2)") - call reset_arrays(5000) - xvalues = linspace(1.0d0, 3.0d0, 5000) - ax = 0.0d0 - bx = xvalues - call integrate_ode_rk( & - xvalues, ax, bx, 5000, yinit=2.0d0, yvalues=yactual, new_xvalues=xrkvalues & - ) + xvalues = linspace(1.0d0, 3.0d0, pts) yexpect = 0.5d0 * xvalues**2 + 1.5d0 - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues, xrkvalues, tolerance=TOL) - @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=1.0d-8) - end subroutine test_integration_simple - - - @test - subroutine test_integration_simple_adaptive() - call set_name("integration (dy=x, x=[1, 3], y0=2, adaptive)") - call reset_arrays(5000) - xvalues = linspace(1.0d0, 3.0d0, 5000) - ax = 0.0d0 - bx = xvalues - call integrate_ode_rk( & - xvalues, & - ax, & - bx, & - 5000, & + call integrate_ode_rk45( & + x0=1.0d0, & + x1=3.0d0, & + ax_func=ax, & + bx_func=bx, & + nbpoints=pts, & yinit=2.0d0, & - yvalues=yactual, & - adaptive=.true., & - new_xvalues=xrkvalues & + yvalues=yactual & ) - yexpect = 0.5d0 * xrkvalues**2 + 1.5d0 - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues(1), xrkvalues(1), tolerance=TOL) - @assertEqual(xvalues(5000), xrkvalues(5000), tolerance=TOL) + @assertEqual(pts, size(yactual)) @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=1.0d-8) - end subroutine test_integration_simple_adaptive + @assertEqual(yexpect, yactual, tolerance=TOL) + + contains + + real(dp) function ax(x) + real(dp), intent(in) :: x + ax = 0.0d0 + end function ax + + real(dp) function bx(x) + real(dp), intent(in) :: x + bx = x + end function bx + end subroutine test_integration_simple @test subroutine test_integration_ode1() call set_name("integration (dy=x**3-xy, x=[0, 2], y0=0)") - call reset_arrays(5000) - xvalues = linspace(0.0d0, 2.0d0, 5000) - ax = -xvalues - bx = xvalues**3 + xvalues = linspace(0.0d0, 2.0d0, pts) yinit = 0.0d0 - call integrate_ode_rk( & - xvalues, ax, bx, 5000, yinit=0.0d0, yvalues=yactual, new_xvalues=xrkvalues & - ) yexpect = xvalues**2 - 2.0d0 + 2.0d0 * exp(-0.5d0 * xvalues**2) - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues, xrkvalues, tolerance=TOL) - @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=1.0d-7) - end subroutine test_integration_ode1 - - - @test - subroutine test_integration_ode1_adaptive() - call set_name("integration (dy=x**3-xy, x=[0, 2], y0=0, adaptive)") - call reset_arrays(5000) - xvalues = linspace(0.0d0, 2.0d0, 5000) - ax = -xvalues - bx = xvalues**3 - yinit = 0.0d0 - call integrate_ode_rk( & - xvalues, & - ax, & - bx, & - 5000, & + call integrate_ode_rk45( & + x0=0.0d0, & + x1=2.0d0, & + ax_func=ax, & + bx_func=bx, & + nbpoints=pts, & yinit=0.0d0, & yvalues=yactual, & - adaptive=.true., & - new_xvalues=xrkvalues & + xvalues=xactual & ) - yexpect = xrkvalues**2 - 2.0d0 + 2.0d0 * exp(-0.5d0 * xrkvalues**2) - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues(1), xrkvalues(1), tolerance=TOL) - @assertEqual(xvalues(5000), xrkvalues(5000), tolerance=TOL) + @assertEqual(pts, size(yactual)) @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=1.0d-7) - end subroutine test_integration_ode1_adaptive - - - @test - subroutine test_integration_ode1_resample() - call set_name("integration (dy=x**3-xy, x=[0, 2], y0=0, resampling)") - call reset_arrays(5000) - xvalues = linspace(0.0d0, 2.0d0, 5000) - ax = -xvalues - bx = xvalues**3 - yinit = 0.0d0 - ! integration will resample from 5000 to 10000 points - deallocate(yactual, xrkvalues, yexpect) - allocate(yactual(10000), xrkvalues(10000), yexpect(10000)) - call integrate_ode_rk( & - xvalues, ax, bx, 10000, yinit=0.0d0, yvalues=yactual, new_xvalues=xrkvalues & - ) - yexpect = xrkvalues**2 - 2.0d0 + 2.0d0 * exp(-0.5d0 * xrkvalues**2) - @assertEqual(5000, size(xvalues)) - @assertEqual(xvalues(1), xrkvalues(1), tolerance=TOL) - @assertEqual(xvalues(5000), xrkvalues(10000), tolerance=TOL) - @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=1.0d-7) - end subroutine test_integration_ode1_resample + @assertEqual(yexpect, yactual, tolerance=TOL) + @assertEqual(pts, size(xactual)) + @assertEqual(xvalues(1), xactual(1), tolerance=TOL) + @assertEqual(xvalues, xactual, tolerance=TOL) + + contains + + real(dp) function ax(x) + real(dp), intent(in) :: x + ax = -x + end function ax + + real(dp) function bx(x) + real(dp), intent(in) :: x + bx = x**3 + end function bx + end subroutine test_integration_ode1 @test subroutine test_integration_ode2() call set_name("integration (dy=-cos(x)y+2x*exp(-sin(x), x=[-2, 3], y0=0)") - call reset_arrays(5000) xvalues = linspace(-2.0d0, 3.0d0, 5000) - ax = -cos(xvalues) - bx = 2.0d0 * xvalues * exp(-sin(xvalues)) - call integrate_ode_rk( & - xvalues, ax, bx, 5000, yinit=0.0d0, yvalues=yactual, new_xvalues=xrkvalues & - ) yexpect = (xvalues**2 - 4.0d0) * exp(-sin(xvalues)) - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues, xrkvalues, tolerance=TOL) - @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=5.0d-6) - end subroutine test_integration_ode2 - - - @test - subroutine test_integration_ode2_adaptive() - call set_name("integration (dy=-cos(x)y+2x*exp(-sin(x), x=[-2, 3], y0=0, adaptive)") - call reset_arrays(5000) - xvalues = linspace(-2.0d0, 3.0d0, 5000) - ax = -cos(xvalues) - bx = 2.0d0 * xvalues * exp(-sin(xvalues)) - call integrate_ode_rk( & - xvalues, & - ax, & - bx, & - 5000, & + call integrate_ode_rk45( & + x0=-2.0d0, & + x1=3.0d0, & + ax_func=ax, & + bx_func=bx, & + nbpoints=pts, & yinit=0.0d0, & - yvalues=yactual, & - adaptive=.true., & - new_xvalues=xrkvalues & + yvalues=yactual & ) - yexpect = (xvalues**2 - 4.0d0) * exp(-sin(xvalues)) - @assertEqual(size(xvalues), size(yactual)) - @assertEqual(xvalues(1), xrkvalues(1), tolerance=TOL) - @assertEqual(xvalues(5000), xrkvalues(5000), tolerance=TOL) + @assertEqual(pts, size(yactual)) @assertEqual(yexpect(1), yactual(1), tolerance=TOL) - @assertEqual(yexpect, yactual, tolerance=5.0d-6) - end subroutine test_integration_ode2_adaptive + @assertEqual(yexpect, yactual, tolerance=TOL) + + contains + + real(dp) function ax(x) + real(dp), intent(in) :: x + ax = -cos(x) + end function ax + + real(dp) function bx(x) + real(dp), intent(in) :: x + bx = 2.0_dp * x * exp(-sin(x)) + end function bx + end subroutine test_integration_ode2 end module mod_test_integration