From b64377977cc3e0aa38acbd49cb9095d5d34194a0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 13 May 2021 16:11:32 -0400 Subject: [PATCH] Started building an objective function for the KE constraint using the lambda function class. --- Makefile | 2 +- src/modules/lambda_function.f90 | 4 ++-- src/modules/module_symba.f90 | 16 ++++++++++++++++ src/symba/symba_frag_pos.f90 | 6 +++--- 4 files changed, 22 insertions(+), 6 deletions(-) diff --git a/Makefile b/Makefile index 19be45a80..4be91714d 100644 --- a/Makefile +++ b/Makefile @@ -47,6 +47,7 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ user.f90 \ swiftest_data_structures.f90 \ + lambda_function.f90\ module_swifter.f90 \ module_helio.f90 \ module_nrutil.f90 \ @@ -54,7 +55,6 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ module_swiftestalloc.f90 \ io.f90 \ module_interfaces.f90 \ - lambda_function.f90\ swiftest.f90 include Makefile.Defines diff --git a/src/modules/lambda_function.f90 b/src/modules/lambda_function.f90 index fcd8c2783..fa672a1f3 100644 --- a/src/modules/lambda_function.f90 +++ b/src/modules/lambda_function.f90 @@ -13,11 +13,11 @@ module lambda_function end type abstract interface - function lambda0(x) + function lambda0(x) result(y) ! Template for a 0 argument function import DP real(DP), intent(in) :: x - real(DP) :: lambda + real(DP) :: y end function end interface diff --git a/src/modules/module_symba.f90 b/src/modules/module_symba.f90 index db506e1ad..fb77a34fb 100644 --- a/src/modules/module_symba.f90 +++ b/src/modules/module_symba.f90 @@ -27,6 +27,7 @@ MODULE module_symba USE swiftest_globals USE module_helio + use lambda_function IMPLICIT NONE INTEGER(I4B), PARAMETER :: NENMAX = 32767 @@ -101,6 +102,21 @@ MODULE module_symba type(swiftest_particle_info), dimension(:), allocatable :: info end type symba_merger + + type, public, extends(lambda_obj) :: symba_vel_lambda_obj + procedure(abstract_objective_func), pointer, nopass :: ke_objective_func_ptr => null() + end type symba_vel_lambda_obj + + abstract interface + function abstract_objective_func(v, m, rhat, t, G, B, L) result(fnorm) + ! Template for the kinetic energy constraint function used for minimizing + import DP + real(DP), dimension(:), intent(in) :: v, m, t, G + real(DP), dimension(:,:), intent(in) :: rhat + real(DP), intent(in) :: B, L + real(DP) :: fnorm + end function + end interface END MODULE module_symba !********************************************************************************************************************************** ! diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index b22a112c9..396a0bf5d 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -370,7 +370,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr end do if (Lambda > 0.0_DP) then lmerge = .false. - v_r_mag(:) = symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambda, tau) + v_r_mag(:) = symba_frag_pos_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau) do i = 1, nfrag v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i) end do @@ -385,12 +385,12 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr return end subroutine symba_frag_pos_kinetic_energy - function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambda, tau) result(v_r_mag) + + function symba_frag_pos_fragment_velocity(nfrag, m_frag, v_r_unit, Lambda, tau) result(v_r_mag) implicit none ! Arguments integer(I4B), intent(in) :: nfrag real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(in) :: x_frag !! Fragment position in the center of mass frame real(DP), dimension(:,:), intent(in) :: v_r_unit !! Radial velocity unit vector for each fragment real(DP), intent(in) :: Lambda !! Sum of the radial kinetic energy of all fragments real(DP), dimension(:), intent(in) :: tau !! Sum of the tangential momentum vector of all fragments