diff --git a/Makefile.Defines b/Makefile.Defines index 291f2c604..07126f842 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -71,7 +71,7 @@ FORTRAN = ifort #AR = xiar #FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) $(GMEM) +#FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 688ac1a6e..d879eec53 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -314,6 +314,7 @@ subroutine calculate_system_energy(linclude_fragments) logical, dimension(:), allocatable :: ltmp logical :: lk_plpl class(swiftest_nbody_system), allocatable :: tmpsys + class(swiftest_parameters), allocatable :: tmpparam ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by @@ -334,27 +335,16 @@ subroutine calculate_system_energy(linclude_fragments) npl_new = npl end if call setup_construct_system(tmpsys, param) + call tmpsys%tp%setup(0, param) deallocate(tmpsys%cb) allocate(tmpsys%cb, source=cb) + allocate(tmpparam, source=param) allocate(ltmp(npl_new)) ltmp(:) = .false. ltmp(1:npl) = .true. call tmpsys%pl%setup(npl_new, param) call tmpsys%pl%fill(pl, ltmp) - tmpsys%pl%mass(1:npl) = tmpsys%pl%mass(1:npl) / mscale - tmpsys%pl%Gmass(1:npl) = tmpsys%pl%Gmass(1:npl) * tscale**2 / rscale**3 - tmpsys%cb%mass = tmpsys%cb%mass / mscale - tmpsys%cb%Gmass = tmpsys%cb%Gmass * tscale**2 / rscale**3 - tmpsys%cb%radius = tmpsys%cb%radius / rscale - tmpsys%pl%radius(1:npl) = tmpsys%pl%radius(1:npl) / rscale - tmpsys%pl%xh(:,1:npl) = tmpsys%pl%xh(:,1:npl) / rscale - tmpsys%pl%vh(:,1:npl) = tmpsys%pl%vh(:,1:npl) / vscale - tmpsys%pl%xb(:,1:npl) = tmpsys%pl%xb(:,1:npl) / rscale - tmpsys%pl%vb(:,1:npl) = tmpsys%pl%vb(:,1:npl) / vscale - tmpsys%pl%rot(:,1:npl) = tmpsys%pl%rot(:,1:npl) * tscale - tmpsys%cb%xb(:) = tmpsys%cb%xb(:) / rscale - tmpsys%cb%vb(:) = tmpsys%cb%vb(:) / vscale - + call tmpsys%rescale(tmpparam, mscale, rscale, tscale) if (linclude_fragments) then ! Append the fragments if they are included ! Energy calculation requires the fragments to be in the system barcyentric frame, s diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 7f4fd140e..4d0e98704 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -306,6 +306,7 @@ module swiftest_classes procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => util_set_msys !! Sets the value of msys from the masses of system bodies. procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum + procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value end type swiftest_nbody_system @@ -984,6 +985,13 @@ end subroutine util_fill_arr_logical end interface interface + module subroutine util_rescale_system(self, param, mscale, dscale, tscale) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters. Returns with new values of the scale vactors and GU + real(DP), intent(in) :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. + end subroutine util_rescale_system + module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1) use lambda_function implicit none @@ -1003,38 +1011,6 @@ module subroutine util_peri_tp(self, system, param) end subroutine util_peri_tp end interface - interface util_solve_linear_system - module function util_solve_linear_system_d(A,b,n,lerr) result(x) - implicit none - integer(I4B), intent(in) :: n - real(DP), dimension(:,:), intent(in) :: A - real(DP), dimension(:), intent(in) :: b - logical, intent(out) :: lerr - real(DP), dimension(n) :: x - end function util_solve_linear_system_d - - module function util_solve_linear_system_q(A,b,n,lerr) result(x) - implicit none - integer(I4B), intent(in) :: n - real(QP), dimension(:,:), intent(in) :: A - real(QP), dimension(:), intent(in) :: b - logical, intent(out) :: lerr - real(QP), dimension(n) :: x - end function util_solve_linear_system_q - end interface - - interface - module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) - use lambda_function - implicit none - class(lambda_obj), intent(inout) :: f !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set - real(DP), dimension(:), intent(in) :: y0in !! Initial value at t=0 - real(DP), intent(in) :: t1 !! Final time - real(DP), intent(in) :: dt0 !! Initial step size guess - real(DP), intent(in) :: tol !! Tolerance on solution - real(DP), dimension(:), allocatable :: y1 !! Final result - end function util_solve_rkf45 - end interface interface util_resize module subroutine util_resize_arr_char_string(arr, nnew) @@ -1142,6 +1118,39 @@ module subroutine util_set_rhill_approximate(self,cb) end subroutine util_set_rhill_approximate end interface + interface util_solve_linear_system + module function util_solve_linear_system_d(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(DP), dimension(:,:), intent(in) :: A + real(DP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(DP), dimension(n) :: x + end function util_solve_linear_system_d + + module function util_solve_linear_system_q(A,b,n,lerr) result(x) + implicit none + integer(I4B), intent(in) :: n + real(QP), dimension(:,:), intent(in) :: A + real(QP), dimension(:), intent(in) :: b + logical, intent(out) :: lerr + real(QP), dimension(n) :: x + end function util_solve_linear_system_q + end interface + + interface + module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1) + use lambda_function + implicit none + class(lambda_obj), intent(inout) :: f !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set + real(DP), dimension(:), intent(in) :: y0in !! Initial value at t=0 + real(DP), intent(in) :: t1 !! Final time + real(DP), intent(in) :: dt0 !! Initial step size guess + real(DP), intent(in) :: tol !! Tolerance on solution + real(DP), dimension(:), allocatable :: y1 !! Final result + end function util_solve_rkf45 + end interface + interface util_sort module subroutine util_sort_i4b(arr) implicit none diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 new file mode 100644 index 000000000..061ecf9a5 --- /dev/null +++ b/src/util/util_rescale.f90 @@ -0,0 +1,51 @@ +submodule (swiftest_classes) s_util_rescale + use swiftest +contains + module subroutine util_rescale_system(self, param, mscale, dscale, tscale) + !! author: David A. Minton + !! + !! Rescales an nbody system to a new set of units. Inputs are the multipliers on the mass (mscale), distance (dscale), and time units (tscale). + !! Rescales all united quantities in the system, as well as the mass conversion factors, gravitational constant, and Einstein's constant in the parameter object. + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters. Returns with new values of the scale vactors and GU + real(DP), intent(in) :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. + ! Internals + real(DP) :: vscale + + param%MU2KG = param%MU2KG * mscale + param%DU2M = param%DU2M * dscale + param%TU2S = param%TU2S * tscale + + ! Calculate the G for the system units + param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) + + ! Calculate the inverse speed of light in the system units + param%inv_c2 = einsteinC * param%TU2S / param%DU2M + param%inv_c2 = (param%inv_c2)**(-2) + + vscale = dscale / tscale + + associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) + + cb%mass = cb%mass / mscale + cb%Gmass = param%GU * cb%mass + cb%radius = cb%radius / dscale + cb%xb(:) = cb%xb(:) / dscale + cb%vb(:) = cb%vb(:) / vscale + pl%mass(1:npl) = pl%mass(1:npl) / mscale + pl%Gmass(1:npl) = param%GU * pl%mass(1:npl) + pl%radius(1:npl) = pl%radius(1:npl) / dscale + pl%xh(:,1:npl) = pl%xh(:,1:npl) / dscale + pl%vh(:,1:npl) = pl%vh(:,1:npl) / vscale + pl%xb(:,1:npl) = pl%xb(:,1:npl) / dscale + pl%vb(:,1:npl) = pl%vb(:,1:npl) / vscale + pl%rot(:,1:npl) = pl%rot(:,1:npl) * tscale + + end associate + + + return + end subroutine util_rescale_system + +end submodule s_util_rescale \ No newline at end of file diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index c6d5aa34f..80d87209c 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -1,6 +1,7 @@ submodule (swiftest_classes) s_util_resize use swiftest contains + module subroutine util_resize_arr_char_string(arr, nnew) !! author: David A. Minton !!