Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Added lambda function class from the Fragmentation branch to use here…
Browse files Browse the repository at this point in the history 
… for the rkf solver
  • Loading branch information
@daminton
daminton committed Jul 13, 2021
1 parent 948670b commit e5ec051
Show file tree
Hide file tree
Showing 4 changed files with 329 additions and 0 deletions.
1 change: 1 addition & 0 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ SWIFTEST_MODULES = swiftest_globals.f90 \
helio_classes.f90 \
symba_classes.f90 \
module_nrutil.f90 \
lambda_function.f90\
swiftest.f90


Expand Down
233 changes: 233 additions & 0 deletions src/modules/lambda_function.f90
View file
Original file line number Diff line number Diff line change
@@ -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

1 change: 1 addition & 0 deletions src/modules/swiftest.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
94 changes: 94 additions & 0 deletions src/util/util_solve.f90
View file
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit e5ec051

Please sign in to comment.