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

Commit

Permalink
Fixed integrator identification problem in the parameter file reader
Browse files Browse the repository at this point in the history
  • Loading branch information
daminton committed Jul 12, 2021
1 parent aa4a21a commit da929af
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 37 deletions.
32 changes: 16 additions & 16 deletions src/gr/gr.f90
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ module pure subroutine gr_getaccb_ns_body(self, system, param)
real(DP) :: rmag, rdotv, vmag2
integer(I4B) :: i

associate(n => self%nbody, cb => system%cb, inv_c2 => param%inv_c2)
associate(n => self%nbody, cb => system%cb, inv_c2 => param%inv_c2)
if (n == 0) return
do i = 1, n
rmag = norm2(self%xh(:,i))
Expand Down Expand Up @@ -92,10 +92,10 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
! Internals
real(DP) :: vmag2, rmag, grterm

associate(c2 => param%inv_c2)
associate(inv_c2 => param%inv_c2)
vmag2 = dot_product(pv(:), pv(:))
rmag = norm2(xh(:))
grterm = 1.0_DP - c2 * (0.5_DP * vmag2 + 3 * mu / rmag)
grterm = 1.0_DP - inv_c2 * (0.5_DP * vmag2 + 3 * mu / rmag)
vh(:) = pv(:) * grterm
end associate
return
Expand Down Expand Up @@ -128,52 +128,52 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
integer(I4B), parameter :: MAXITER = 50
real(DP),parameter :: TOL = 1.0e-12_DP

associate (c2 => param%inv_c2)
associate(inv_c2 => param%inv_c2)
pv(1:NDIM) = vh(1:NDIM) ! Initial guess
rterm = 3 * mu / norm2(xh(:))
v2 = dot_product(vh(:), vh(:))
do n = 1, MAXITER
pv2 = dot_product(pv(:), pv(:))
G = 1.0_DP - c2 * (0.5_DP * pv2 + rterm)
G = 1.0_DP - inv_c2 * (0.5_DP * pv2 + rterm)
F(:) = pv(:) * G - vh(:)
if (abs(sum(F) / v2 ) < TOL) exit ! Root found

! Calculate the Jacobian
do k = 1, NDIM
do i = 1, NDIM
if (i == k) then
J(i,k) = G - c2 * pv(k)
J(i,k) = G - inv_c2 * pv(k)
else
J(i,k) = - c2 * pv(k)
J(i,k) = -inv_c2 * pv(k)
end if
end do
end do

! Inverse of the Jacobian
det = J(1,1) * (J(3,3) * J(2,2) - J(3,2) * J(2,3))
det = det - J(2,1) * (J(3,3) * J(1,2)-J(3,2) * J(1,3))
det = det + J(3,1) * (J(2,3) * J(1,2)-J(2,2) * J(1,3))

Jinv(1,1) = J(3,3) * J(2,2) - J(3,2) * J(2,3)
Jinv(1,2) = -(J(3,3) * J(1,2) - J(3,2) * J(1,3))
Jinv(1,3) = J(2,3) * J(1,2) - J(2,2) * J(1,3)

Jinv(2,1) = -(J(3,3) * J(2,1) - J(3,1) * J(2,3))
Jinv(2,2) = J(3,3) * J(1,1) - J(3,1) * J(1,3)
Jinv(2,3) = -(J(2,3) * J(1,1) - J(2,1) * J(1,3))

Jinv(3,1) = J(3,2) * J(2,1) - J(3,1) * J(2,2)
Jinv(3,2) = -(J(3,2) * J(1,1) - J(3,1) * J(1,2))
Jinv(3,3) = J(2,2) * J(1,1) - J(2,1) * J(1,2)

Jinv = Jinv * det

do i = 1, NDIM
pv(i) = pv(i) - dot_product(Jinv(i,:) ,F(:))
pv(i) = pv(i) - dot_product(Jinv(i,:), F(:))
end do
end do

end associate

return
end subroutine gr_vel2pseudovel

Expand Down
39 changes: 20 additions & 19 deletions src/io/io.f90
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
character(STRMAX) :: line !! Line of the input file
character (len=:), allocatable :: line_trim,param_name, param_value !! Strings used to parse the param file
character(*),parameter :: linefmt = '(A)' !! Format code for simple text string
integer(I4B) :: integrator !! Symbolic name of integrator being used

integrator = v_list(1)
! Parse the file line by line, extracting tokens then matching them up with known parameters if possible

do
read(unit = unit, fmt = linefmt, iostat = iostat, end = 1) line
line_trim = trim(adjustl(line))
Expand Down Expand Up @@ -214,7 +213,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
write(*,*) "ENC_OUT = ",trim(adjustl(self%encounter_file))
write(*,*) "EXTRA_FORCE = ",self%lextra_force
write(*,*) "BIG_DISCARD = ",self%lbig_discard
if (self%lenergy) write(*,*) "ENERGY = ",self%lenergy
write(*,*) "ENERGY = ",self%lenergy
write(*,*) "MU2KG = ",self%MU2KG
write(*,*) "TU2S = ",self%TU2S
write(*,*) "DU2M = ",self%DU2M
Expand All @@ -232,21 +231,23 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
self%inv_c2 = einsteinC * self%TU2S / self%DU2M
self%inv_c2 = (self%inv_c2)**(-2)

if (integrator == RMVS) then
if (.not.self%lclose) then
write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.'
iostat = -1
return
associate(integrator => v_list(1))
if (integrator == RMVS) then
if (.not.self%lclose) then
write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.'
iostat = -1
return
end if
end if
end if

! Determine if the GR flag is set correctly for this integrator
select case(integrator)
case(WHM, RMVS)
write(*,*) "GR = ", self%lgr
case default
write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.'
end select

! Determine if the GR flag is set correctly for this integrator
select case(integrator)
case(WHM, RMVS)
write(*,*) "GR = ", self%lgr
case default
write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.'
end select
end associate

iostat = 0

Expand Down Expand Up @@ -682,8 +683,8 @@ module subroutine io_read_param_in(self, param_file_name)
!! Adapted from Martin Duncan's Swift routine io_init_param.f
implicit none
! Arguments
class(swiftest_parameters),intent(out) :: self !! Current run configuration parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
class(swiftest_parameters),intent(inout) :: self !! Current run configuration parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
! Internals
integer(I4B), parameter :: LUN = 7 !! Unit number of input file
integer(I4B) :: ierr = 0 !! Input error code
Expand Down
4 changes: 2 additions & 2 deletions src/modules/swiftest_classes.f90
Original file line number Diff line number Diff line change
Expand Up @@ -488,8 +488,8 @@ end subroutine io_read_cb_in

module subroutine io_read_param_in(self, param_file_name)
implicit none
class(swiftest_parameters), intent(out) :: self !! Current run configuration parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
class(swiftest_parameters), intent(inout) :: self !! Current run configuration parameters
character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
end subroutine io_read_param_in

module subroutine io_read_frame_body(self, iu, param, form, ierr)
Expand Down

0 comments on commit da929af

Please sign in to comment.