diff --git a/Makefile b/Makefile index e78f0ceeb..7154c63b1 100644 --- a/Makefile +++ b/Makefile @@ -52,6 +52,7 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ helio_classes.f90 \ symba_classes.f90 \ module_nrutil.f90 \ + lambda_function.f90\ swiftest.f90 diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90 new file mode 100644 index 000000000..febfd457d --- /dev/null +++ b/src/modules/lambda_function.f90 @@ -0,0 +1,233 @@ +module lambda_function + !! author: David A. Minton + !! + !! Defines a class that can enable objects that behave like lambda functions. + !! + !! To use this class, define a type of either lambda_obj or lambda_obj_err, or extend the lambda_obj class as necessary, such that an interface that matches the function you wish to lambdafy. + !! Once defined, the lambda object can evaluate itself by calling the type-bound procedure eval. e.g. f%eval(x) (or f%eval(x, lerr), f%eval(x, [argument list], etc)) + !! + !! ******************************************************************************************************************************************************************************************** + !! Example - Defining a lambda function f(x,rval,ival) where rval and ival are a real and integer argument, respectively. This implementation uses an abstract interface, though this is not + !! strictly necessary unless you want to bind more than one function with the same interface. + !! ******************************************************************************************************************************************************************************************** + !! + !! module lambda_new + !! use swiftest ! This will bring in the lambda_function module + !! ! Define types in a module + !! + !! type, extends(lambda_obj) :: lambda_obj_ri_args + !! procedure(abstract_lambda_ri_args), pointer, nopass :: lambdaptr_ri_args => null() + !! real(DP) :: rval !! Real parameter + !! integer(I4B) :: ival !! Integer paramete + !! contains + !! generic :: init => lambda_ri_args_init + !! procedure :: eval => lambda_ri_args_eval + !! procedure, nopass :: lambda_ri_args_init + !! final :: lambda_ri_args_destroy + !! end type + !! interface lambda_obj + !! module procedure lambda_ri_args_init + !! end interface + !! + !! abstract interface + !! function abstract_lambda_ri_args(x, rval, ival) result(y) + !! !Template for the lambda function + !! import DP, I4B + !! real(DP), dimension(:), intent(in) :: x !! Dependent variable + !! real(DP), intent(in) :: rval !! Real parameter + !! integer(I4B), intent(in) :: ival !! Integer parameter + !! real(DP) :: y !! Real result + !! end function + !! end interface + !! + !! contains + !! type(lambda_obj_ri_args) function lambda_ri_args_init(lambda, rval, ival) + !! !! Initializes the lambda function parameters (can be used as a structure constructor) + !! implicit none + !! procedure(abstract_lambda_ri_args) :: lambda !! The lambda function that will be passed + !! real(DP), intent(in) :: rval !! Real parameter + !! integer(I4B), intent(in) :: ival !! Integer parameter + !! + !! ! Assign the procedure passed to this function to the procedure pointer + !! lambda_ri_args_init%lambdaptr_ri_args => lambda + !! + !! ! Assign the argument values + !! lambda_ri_args_init%rval = rval + !! lambda_ri_args_init%ival = ival + !! return + !! end function lambda_ri_args_init + !! + !! function lambda_ri_args_eval(self, x) result(y) + !! !! Defines the evaluation method, allowing the lambda function to be called with a single argument + !! implicit none + !! class(lambda_obj_ri_args), intent(inout) :: self + !! real(DP), dimension(:), intent(in) :: x + !! real(DP) :: y + !! + !! if (associated(self%lambdaptr_ri_args)) then + !! y = self%lambdaptr_ri_args(x, self%rval, self%ival) + !! self%lastval = y + !! if (allocated(self%lastarg)) deallocate(self%lastarg) + !! allocate(self%lastarg, source=x) + !! else + !! error stop "Lambda function was not initialized" + !! end if + !! end function lambda_ri_args_eval + !! + !! subroutine lambda_ri_args_destroy(self) + !! !! Finalizer method. Use this as a template for cleaning up the object upon destruction, such as nullifying pointers + !! implicit none + !! type(lambda_obj_ri_args) :: self + !! if (associated(self%lambdaptr_ri_args)) nullify(self%lambdaptr_ri_args) + !! end subroutine lambda_ri_args_destroy + !! + !! function example_function(x, rval, ival) result(y) + !! !This is the actual function you are going to use as the lambda function. Its interface must match the abstract interface previously defined + !! implicit none + !! ! Arguments + !! real(DP), dimension(:), intent(in) :: x + !! real(DP), intent(in) :: rval + !! integer(I4B), intent(in) :: ival + !! ! Result + !! real(DP) :: y + !! ! Internals + !! integer(I4B) :: i, n + !! n = size(x) + !! y = 42._DP * ival + !! do i = 1, n + !! y = y + x(i)**2 + !! end do + !! return + !! end function example_function + !! end module lambda_new + !! + !! program usage + !! use swiftest + !! use lambda_new + !! implicit none + !! type(lambda_obj_ri_args) :: f + !! real(DP) :: sigma_par + !! integer(I4B) :: iwonky, i,j + !! real(DP), dimension(12) :: xarr + !! + !! sigma_par = 3.14_DP + !! iwonky = 13 + !! + !! f = lambda_obj(example_function, sigma_par, iwonky) + !! do i = 1, 10 + !! xarr(:) = [(j * 0.25_DP / i, j=1, 12)] + !! write(*,*) i,f%eval(xarr) + !! end do + !! end program usage + !! ******************************************************************************************************************************************************************************************** + + use swiftest_globals + implicit none + public + + type :: lambda_obj + !! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers + procedure(lambda0), pointer, nopass :: lambdaptr => null() + real(DP) :: lastval + real(DP),dimension(:), allocatable :: lastarg + contains + generic :: init => lambda_init_0 + procedure :: eval => lambda_eval_0 + procedure, nopass :: lambda_init_0 + final :: lambda_destroy + end type + + type, extends(lambda_obj) :: lambda_obj_err + !! Extended class for an lambda function object. This object takes allows for the return of a logical error flag during evaluation of the function. + procedure(lambda0err), pointer, nopass :: lambdaptr_err => null() + logical :: lerr + contains + generic :: init => lambda_init_0_err + procedure :: eval => lambda_eval_0_err + procedure, nopass :: lambda_init_0_err + end type + interface lambda_obj + module procedure lambda_init_0 + module procedure lambda_init_0_err + end interface + + abstract interface + function lambda0(x) result(y) + ! Template for a 0 argument function + import DP + real(DP), dimension(:), intent(in) :: x + real(DP) :: y + end function + function lambda0err(x, lerr) result(y) + ! Template for a 0 argument function that returns an error value + import DP + real(DP), dimension(:), intent(in) :: x + logical, intent(out) :: lerr + real(DP) :: y + end function + end interface + + contains + type(lambda_obj) function lambda_init_0(lambda) + implicit none + ! Arguments + procedure(lambda0) :: lambda + + lambda_init_0%lambdaptr => lambda + return + end function lambda_init_0 + + type(lambda_obj_err) function lambda_init_0_err(lambda, lerr) + implicit none + ! Arguments + procedure(lambda0err) :: lambda + logical, intent(in) :: lerr + lambda_init_0_err%lambdaptr_err => lambda + lambda_init_0_err%lerr = lerr + return + end function lambda_init_0_err + + function lambda_eval_0(self, x) result(y) + implicit none + ! Arguments + class(lambda_obj), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + ! Result + real(DP) :: y + + if (associated(self%lambdaptr)) then + y = self%lambdaptr(x) + self%lastval = y + if (allocated(self%lastarg)) deallocate(self%lastarg) + allocate(self%lastarg, source=x) + else + error stop "Lambda function was not initialized" + end if + end function lambda_eval_0 + + function lambda_eval_0_err(self, x) result(y) + implicit none + ! Arguments + class(lambda_obj_err), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + ! Result + real(DP) :: y + + if (associated(self%lambdaptr_err)) then + y = self%lambdaptr_err(x, self%lerr) + self%lastval = y + if (allocated(self%lastarg)) deallocate(self%lastarg) + allocate(self%lastarg, source=x) + else + error stop "Lambda function was not initialized" + end if + end function lambda_eval_0_err + + subroutine lambda_destroy(self) + implicit none + type(lambda_obj) :: self + if (associated(self%lambdaptr)) nullify(self%lambdaptr) + end subroutine lambda_destroy + +end module lambda_function + diff --git a/src/modules/swiftest.f90 b/src/modules/swiftest.f90 index fd8ffaf1b..2ab44a9d5 100644 --- a/src/modules/swiftest.f90 +++ b/src/modules/swiftest.f90 @@ -11,6 +11,7 @@ module swiftest use helio_classes use symba_classes use module_nrutil + use lambda_function !use advisor_annotate !$ use omp_lib implicit none diff --git a/src/util/util_solve.f90 b/src/util/util_solve.f90 new file mode 100644 index 000000000..616ef9a4a --- /dev/null +++ b/src/util/util_solve.f90 @@ -0,0 +1,94 @@ +! submodule(swiftest_classes) s_util_solve +! use swiftest +! contains +! subroutine util_solve_rkf45(xv, dt0, tol) +! !! author: David A. Minton +! !! +! !! Implements the 4th order Runge-Kutta-Fehlberg ODE solver for initial value problemswith 5th order adaptive step size control. +! implicit none +! ! Arguments +! real(DP), dimension(:), intent(in) :: xv !! The dependent variable +! real(DP), intent(in) :: dt0, tol ! Output cadence time step (also used as initial step size guess) and error tolerance +! integer :: i, n, nsteps ! The number of steps to generate output +! integer, parameter :: maxredux = 1000 ! Maximum number of times step size can be reduced +! real(DP),dimension(:), allocatable :: y,y0,ynorm ! Internal temporary variable used to store intermediate results until total number of steps is known +! integer, parameter :: rks = 6 ! Number of RK stages +! real(DP),dimension(rks, rks - 1),parameter :: rkf45_btab = reshape( & ! Butcher tableau for Runge-Kutta-Fehlberg method +! (/ 1./4., 1./4., 0., 0., 0., 0.,& +! 3./8., 3./32., 9./32., 0., 0., 0.,& +! 12./13., 1932./2197., -7200./2197., 7296./2197., 0., 0.,& +! 1., 439./216., -8., 3680./513., -845./4104., 0.,& +! 1./2., -8./27., 2., -3544./2565., 1859./4104., -11./40./), shape(rkf45_btab)) +! real(DP),dimension(rks),parameter :: rkf4_coeff = (/ 25./216., 0., 1408./2565. , 2197./4104. , -1./5., 0. /) +! real(DP),dimension(rks),parameter :: rkf5_coeff = (/ 16./135., 0., 6656./12825., 28561./56430., -9./50., 2./55. /) +! real(DP), dimension(:, :), allocatable :: k ! Runge-Kutta coefficient vector +! integer :: rkn ! Runge-Kutta loop index +! integer :: ndim ! Number of dimensions of the problem +! real(DP) :: dt, trem ! Current step size and total time remaining +! real(DP) :: s, yerr, yscale ! Step size reduction factor, error in dependent variable, and error scale factor +! real(DP), parameter :: dtfac = 0.95_DP ! Step size reduction safety factor (Value just under 1.0 to prevent adaptive step size control from discarding steps too aggressively) +! real(DP) :: dtmean ! Mean step size +! integer :: ntot ! Total number of steps (used in mean step size calculation) +! real(DP) :: xscale, vscale + +! ndim = size(xv, 1) +! nsteps = size(xv, 2) +! allocate(k(ndim, rks)) +! allocate(y(ndim)) +! allocate(y0(ndim)) +! allocate(ynorm(ndim)) + +! dt = dt0 +! dtmean = 0.0_DP +! ntot = 0 + +! do n = 2, nsteps +! y0(:) = xv(:, n - 1) +! trem = dt0 +! do +! yscale = norm2(y0(:)) +! xscale = norm2(y0(1:2)) +! vscale = norm2(y0(3:4)) +! do i = 1, maxredux +! do rkn = 1, rks +! y(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1)) +! k(:, rkn) = dt * derivs(y(:)) +! end do +! ! Now determine if the step size needs adjusting +! ynorm(:) = matmul(k(:,:), (rkf5_coeff(:) - rkf4_coeff(:))) +! ynorm(1:2) = ynorm(1:2) / xscale +! ynorm(3:4) = ynorm(3:4) / vscale +! !ynorm(:) = ynorm(:) / yscale +! yerr = norm2(ynorm(:)) +! s = (tol / (2 * yerr))**(0.25_DP) +! dt = min(s * dtfac * dt, trem) ! Alter step size either up or down +! if (s >= 1.0_DP) exit ! Good step! +! if (i == maxredux) then +! write(*,*) 'Something has gone wrong!!' +! stop +! end if +! end do + +! ! Compute new value +! y(:) = y0(:) + matmul(k(:, :), rkf4_coeff(:)) +! trem = trem - dt +! ntot = ntot + 1 +! dtmean = dtmean + dt +! if (trem <= 0._DP) exit +! y0(:) = y(:) +! end do + +! xv(:,n) = y(:) +! end do + +! dtmean = dtmean / ntot +! write(*,*) 'Total number of steps taken: ',ntot +! write(*,*) 'Mean step size: ', dtmean / (2 * pi) + +! deallocate(k,y,y0) + +! return + +! end subroutine util_solve_rkf45 + +! end submodule s_util_solve \ No newline at end of file