From d0126da35c2da980e0ceaf80be926e803257cdeb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Jul 2021 17:23:21 -0400 Subject: [PATCH] Created the structure needed to send spin equations into the adaptive RKF45 solver --- .../swiftest_rmvs_vs_swifter_rmvs.ipynb | 4 +- src/modules/lambda_function.f90 | 46 +++++++- src/tides/tides_spin_step.f90 | 110 ++++++++++++++++++ src/util/util_solve.f90 | 27 ++++- 4 files changed, 176 insertions(+), 11 deletions(-) diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb index 973a1d8c1..d0d223ce7 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb @@ -622,7 +622,7 @@ " 1.17040422e-11])\n", "Coordinates:\n", " id int64 2\n", - " * time (d) (time (d)) float64 0.0 11.0 22.0 ... 3.63e+03 3.641e+03 3.652e+03
    • id
      ()
      int64
      2
      array(2)
    • time (d)
      (time (d))
      float64
      0.0 11.0 ... 3.641e+03 3.652e+03
      array([   0.,   11.,   22., ..., 3630., 3641., 3652.])
    • " ], "text/plain": [ "\n", diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90 index febfd457d..cd7e180cb 100644 --- a/src/modules/lambda_function.f90 +++ b/src/modules/lambda_function.f90 @@ -146,9 +146,19 @@ module lambda_function procedure :: eval => lambda_eval_0_err procedure, nopass :: lambda_init_0_err end type + + type, extends(lambda_obj) :: lambda_obj_tvar + !! 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(lambda0tvar), pointer, nopass :: lambdaptr_tvar => null() + contains + generic :: init => lambda_init_tvar + procedure :: evalt => lambda_eval_tvar + procedure, nopass :: lambda_init_tvar + end type interface lambda_obj module procedure lambda_init_0 module procedure lambda_init_0_err + module procedure lambda_init_tvar end interface abstract interface @@ -158,6 +168,7 @@ function lambda0(x) result(y) 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 @@ -165,6 +176,14 @@ function lambda0err(x, lerr) result(y) logical, intent(out) :: lerr real(DP) :: y end function + + function lambda0tvar(x, t) result(y) + ! Template for a 0 argument function that returns an error value + import DP + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + real(DP), dimension(:), allocatable :: y + end function end interface contains @@ -172,7 +191,6 @@ 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 @@ -186,6 +204,15 @@ type(lambda_obj_err) function lambda_init_0_err(lambda, lerr) lambda_init_0_err%lerr = lerr return end function lambda_init_0_err + + type(lambda_obj_tvar) function lambda_init_tvar(lambda, t) + implicit none + ! Arguments + procedure(lambda0tvar) :: lambda + real(DP), intent(in) :: t + lambda_init_tvar%lambdaptr_tvar => lambda + return + end function lambda_init_tvar function lambda_eval_0(self, x) result(y) implicit none @@ -194,7 +221,6 @@ function lambda_eval_0(self, x) result(y) real(DP), dimension(:), intent(in) :: x ! Result real(DP) :: y - if (associated(self%lambdaptr)) then y = self%lambdaptr(x) self%lastval = y @@ -212,7 +238,6 @@ function lambda_eval_0_err(self, x) result(y) 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 @@ -223,6 +248,21 @@ function lambda_eval_0_err(self, x) result(y) end if end function lambda_eval_0_err + function lambda_eval_tvar(self, x, t) result(y) + implicit none + ! Arguments + class(lambda_obj_tvar), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + ! Result + real(DP), dimension(:), allocatable :: y + if (associated(self%lambdaptr_tvar)) then + y = self%lambdaptr_tvar(x,t) + else + error stop "Lambda function was not initialized" + end if + end function lambda_eval_tvar + subroutine lambda_destroy(self) implicit none type(lambda_obj) :: self diff --git a/src/tides/tides_spin_step.f90 b/src/tides/tides_spin_step.f90 index 52f2bae0f..d6029a7e4 100644 --- a/src/tides/tides_spin_step.f90 +++ b/src/tides/tides_spin_step.f90 @@ -1,15 +1,125 @@ submodule(swiftest_classes) s_tides_step_spin use swiftest + + type, extends(lambda_obj_tvar) :: tides_derivs_func + !! 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(tidederiv), pointer, nopass :: lambdaptr_tides_deriv + real(DP), dimension(:,:), allocatable :: xbeg + real(DP), dimension(:,:), allocatable :: xend + real(DP) :: dt + contains + generic :: init => tides_derivs_init + procedure :: evalt => tides_derivs_eval + procedure, nopass :: tides_derivs_init + end type + interface lambda_obj + module procedure tides_derivs_init + end interface + abstract interface + function tidederiv(x, t, dt, xbeg, xend) result(y) + ! Template for a 0 argument function + import DP, swiftest_nbody_system + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + real(DP), intent(in) :: dt + real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: xend + real(DP), dimension(:), allocatable :: y + end function + end interface + contains module subroutine tides_step_spin_system(self, param, t, dt) !! author: Jennifer L.L. Pouplin and David A. Minton !! !! Integrates the spin equations for central and massive bodies of the system subjected to tides. implicit none + ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Simulation time real(DP), intent(in) :: dt !! Current stepsize + ! Internals + real(DP), dimension(:), allocatable :: rot0, rot1 + real(DP) :: subt + real(DP), parameter :: tol=1e-6_DP !! Just a guess at the moment + real(DP) :: subdt + + associate(pl => self%pl, npl => self%pl%nbody, cb => self%cb) + allocate(rot0(NDIM*(npl+1))) + rot0 = [pack(pl%rot(:,1:npl),.true.), pack(cb%rot(:),.true.)] + ! Use this space call the ode_solver, passing tides_spin_derivs as the function: + subdt = dt / 20._DP + !rot1(:) = util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%xbeg, pl%xend), rot0, dt, subdt tol) + ! Recover with unpack + !pl%rot(:,1:npl) = unpack(rot1... + !cb%rot(:) = unpack(rot1... + end associate + return end subroutine tides_step_spin_system + + function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... + !! author: Jennifer L.L. Pouplin and David A. Minton + !! + !! function used to calculate the derivatives that are fed to the ODE solver + implicit none + ! Arguments + real(DP), dimension(:,:), intent(in) :: rot_pl_cb !! Array of rotations. The last element is the central body, and all others are massive bodies + real(DP), intent(in) :: t !! Current time, which is used to interpolate the massive body positions + real(DP), intent(in) :: dt !! Total step size + real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: xend + ! Internals + real(DP), dimension(:,:), allocatable :: drot + real(DP), dimension(:), allocatable :: flatrot + real(DP), dimension(NDIM) :: N_Tcb, N_Rcb, N_Tpl, N_Rpl, xinterp + real(DP) :: C_cb, C_pl, r_dot_rot_cb, r_dot_rot_pl, rmag + integer(I4B) :: i, n + + + n = size(rot_pl_cb,2) + if (allocated(drot)) deallocate(drot) + allocate(drot, mold=rot_pl_cb) + drot(:,:) = 0.0_DP + do i = 1,n-1 + xinterp(:) = xbeg(:,i) + t / dt * (xend(:,i) - xbeg(:,i)) + ! Calculate Ncb and Npl as a function of xinterp + !drot(:,i) = -Mcb / (Mcb + Mpl(i)) * (N_Tpl + N_Rpl) + !drot(:,n) = drot(:,n) - Mcb / (Mcb + Mpl(i) * (N_Tcb + N_Rcb) + ! + end do + + end function tides_spin_derivs + + function tides_derivs_eval(self, x, t) result(y) + implicit none + ! Arguments + class(tides_derivs_func), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + ! Result + real(DP), dimension(:), allocatable :: y + if (associated(self%lambdaptr_tides_deriv)) then + y = self%lambdaptr_tides_deriv(x, t, self%dt, self%xbeg, self%xend) + else + error stop "Lambda function was not initialized" + end if + end function tides_derivs_eval + + function tides_derivs_init(lambda, dt, xbeg, xend) result(f) + implicit none + ! Arguments + procedure(tidederiv) :: lambda + real(DP), intent(in) :: dt + real(DP), dimension(:,:), intent(in) :: xbeg + real(DP), dimension(:,:), intent(in) :: xend + ! Result + type(tides_derivs_func) :: f + f%lambdaptr_tides_deriv => lambda + f%dt = dt + allocate(f%xbeg, source = xbeg) + allocate(f%xend, source = xend) + return + end function tides_derivs_init end submodule s_tides_step_spin \ No newline at end of file diff --git a/src/util/util_solve.f90 b/src/util/util_solve.f90 index 8b1dd0992..255137f3d 100644 --- a/src/util/util_solve.f90 +++ b/src/util/util_solve.f90 @@ -29,10 +29,10 @@ function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) 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 real(DP), dimension(:), allocatable :: ynorm !! Normalized y value used for adaptive step size control - real(DP), dimension(:), allocatable :: y0 !! Value of y at the beginning of each substep + real(DP), dimension(:), allocatable :: y0 !! Value of y at the beginning of each substep integer(I4B) :: Nvar !! Number of variables in problem integer(I4B) :: rkn !! Runge-Kutta loop index - real(DP) :: dt, trem !! Current step size and total time remaining + real(DP) :: t, x1, dt, trem !! Current time, step size and total time remaining real(DP) :: s, yerr, yscale !! Step size reduction factor, error in dependent variable, and error scale factor integer(I4B) :: i, n @@ -45,13 +45,27 @@ function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) dt = dt0 trem = t1 + t = 0._DP do yscale = norm2(y0(:)) do i = 1, MAXREDUX - do rkn = 1, RKS - y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1)) - k(:, rkn) = dt * f%eval(y1(:)) - end do + select type(f) + class is (lambda_obj_tvar) + do rkn = 1, RKS + y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1)) + if (rkn == 1) then + x1 = t + else + x1 = t + rkf45_btab(1,rkn-1) + end if + k(:, rkn) = dt * f%evalt(y1(:), t) + end do + class is (lambda_obj) + do rkn = 1, RKS + y1(:) = y0(:) + matmul(k(:, 1:rkn - 1), rkf45_btab(2:rkn, rkn - 1)) + k(:, rkn) = dt * f%eval(y1(:)) + end do + end select ! Now determine if the step size needs adjusting ynorm(:) = matmul(k(:,:), (rkf5_coeff(:) - rkf4_coeff(:))) / yscale yerr = norm2(ynorm(:)) @@ -67,6 +81,7 @@ function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) ! Compute new value then step ahead in time y1(:) = y0(:) + matmul(k(:, :), rkf4_coeff(:)) trem = trem - dt + t = t + dt if (trem <= 0._DP) exit y0(:) = y1(:) end do