-
Notifications
You must be signed in to change notification settings - Fork 136
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
mixed precision time_interp #1252
Changes from 33 commits
46d8e44
dde0b88
bc9e4d9
2977578
b3374f6
cca0c7d
2795ed2
99980a2
9075302
f2ddd83
2d4c4d8
3e47f37
eab2cba
c50d55c
a0897a5
8878b54
e19a188
5ec7ac6
712330b
673b901
1389016
109aa21
7fb5628
40712be
94434d9
8c72630
75a6b7f
8775ced
608abb7
3a7d6b2
8be0313
81a9fa9
27f31d0
499865b
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -23,19 +23,24 @@ program test_time_interp | |
use time_manager_mod, only: get_date, set_time, set_date, time_manager_init, set_calendar_type, operator(+) | ||
use time_manager_mod, only: JULIAN, time_type, increment_time, NOLEAP, print_date | ||
use time_interp_mod, only: time_interp_init, time_interp, NONE, YEAR, MONTH, DAY | ||
use time_manager_mod, only: operator(<=), operator(>=), operator(==) | ||
use platform_mod | ||
|
||
implicit none | ||
|
||
integer, parameter :: num_Time=6 | ||
integer, parameter :: num_Time=6, kindl = TI_TEST_KIND_ | ||
type(time_type) :: Time_beg, Time_end, Time(num_Time) | ||
type(time_type), allocatable, dimension(:) :: Timelist | ||
integer :: index1, index2, mo, yr, timelist_len, outunit, ntest, nline | ||
real :: weight | ||
integer :: index1, index2, mo, yr, outunit, ntest, nline | ||
real(TI_TEST_KIND_) :: weight | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It seems like the test is a guide on how to use time_interp. Is there a way to check that the answers are as expected? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I added checks for all but the last, since that one cycles through instead of going to set dates. |
||
real(TI_TEST_KIND_) :: ref_weights(num_Time), ref_weights_leap(num_Time) | ||
real(TI_TEST_KIND_), parameter :: SMALL = 1.0e-7_kindl ! r4 will fail with 8 | ||
real(TI_TEST_KIND_), parameter :: midpoint = 0.483870967741935_kindl | ||
real(TI_TEST_KIND_), parameter :: day_before_leap_day = 0.964285714285714_kindl | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. where did these numbers come from? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just reference weights from the calculation, its a fraction of time passed over a interval. So for this its the midway point of a month (jan 16). I think i commented on it a bit lower down. I might be able to replace these with the actual time_type calculation if thats preferred, i guess i just like seeing what the actual value ends up as. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. could the actual time_type calculation be used and add the numbers as comments? Is this perhaps the reason you're not getting exactly agreeing answers? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. After looking into it i don't think i can reproduce the calculation unless i copy in a whole routine (set_modtime). I think this way has some benefits though, like if the time_type arithmentic overrides return the wrong thing a check with the same operations will just be checking wrong answers against wrong answers. |
||
real(TI_TEST_KIND_), parameter :: day_before_leap_day_with_ly = 0.931034482758621_kindl | ||
|
||
integer :: nmin, nmax | ||
|
||
namelist / test_time_interp_nml / timelist_len | ||
|
||
call fms_init | ||
outunit = stdout() | ||
call set_calendar_type(JULIAN) | ||
|
@@ -50,8 +55,23 @@ program test_time_interp | |
Time(5) = set_date(1,12,16) | ||
Time(6) = Time_end | ||
|
||
ref_weights(1) = 0.0 ! on 'edge' (timeList value) | ||
ref_weights(2) = midpoint ! rough midpoint of a month ie. jan 16 | ||
ref_weights(3) = 0.0 | ||
ref_weights(4) = 0.0 | ||
ref_weights(5) = midpoint | ||
ref_weights(6) = 0.0 | ||
|
||
ref_weights_leap(1) = 0.0 ! on 'edge' (timeList value) | ||
ref_weights_leap(2) = day_before_leap_day ! feb 28th | ||
ref_weights_leap(3) = midpoint | ||
ref_weights_leap(4) = 0.0 | ||
ref_weights_leap(5) = day_before_leap_day | ||
ref_weights_leap(6) = day_before_leap_day ! checks that 29th gives same result | ||
|
||
! Tests with modulo time | ||
do nline=1,3 | ||
|
||
if(nline == 1) then | ||
allocate(Timelist(12)) | ||
do mo=1,12 | ||
|
@@ -72,6 +92,7 @@ program test_time_interp | |
endif | ||
|
||
do ntest=1,num_Time | ||
print *, ntest | ||
call diagram(nline,ntest,modulo_time=.true.) | ||
call time_interp(Time(ntest), Time_beg, Time_end, Timelist, weight, index1, index2) | ||
write(outunit,*) 'time_interp_modulo:' | ||
|
@@ -84,17 +105,30 @@ program test_time_interp | |
write(outunit,99) index1,index2,weight | ||
write(outunit,'()') | ||
|
||
if(.not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) & | ||
call mpp_error(FATAL, "test_time_interp: invalid indices from time_interp_timelist") | ||
if(abs(weight - ref_weights(ntest)) .gt. SMALL) & | ||
call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference") | ||
|
||
call time_interp(Time(ntest), Timelist, weight, index1, index2, modtime=YEAR) | ||
write(outunit,*) 'time_interp_list with modtime=YEAR:' | ||
write(outunit,'()') | ||
call print_date(Time(ntest), 'Time =') | ||
call print_date(Timelist(1), 'Timelist(1)=') | ||
call print_date(Timelist(size(Timelist(:))),'Timelist(n)=') | ||
write(outunit,99) index1,index2,weight | ||
|
||
if(.not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) & | ||
call mpp_error(FATAL, "test_time_interp: invalid indices from time_interp_modulo") | ||
if(abs(weight - ref_weights(ntest)) .gt. SMALL) & | ||
call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference") | ||
|
||
enddo | ||
deallocate(Timelist) | ||
enddo | ||
|
||
|
||
|
||
! Tests without modulo time | ||
do nline=1,3 | ||
if(nline == 1) then | ||
|
@@ -132,6 +166,12 @@ program test_time_interp | |
call print_date(Timelist(1), 'Timelist(1)=') | ||
call print_date(Timelist(size(Timelist(:))),'Timelist(n)=') | ||
write(outunit,99) index1,index2,weight | ||
|
||
if( .not. is_valid_indices(index1, index2, TimeList, Time(ntest), weight, NONE)) & | ||
call mpp_error(FATAL, "invalid result without modtime") | ||
if(abs(weight - ref_weights(ntest)) .gt. SMALL) & | ||
call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference") | ||
|
||
enddo | ||
deallocate(Timelist) | ||
enddo | ||
|
@@ -171,9 +211,20 @@ program test_time_interp | |
call print_date(Time(ntest),' Time =') | ||
write(outunit,99) index1,index2,weight | ||
write(outunit,'()') | ||
if( .not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) & | ||
call mpp_error(FATAL, 'invalid results for indices with leap year correction') | ||
if(abs(weight - ref_weights_leap(ntest)) .gt. SMALL) & | ||
call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference") | ||
enddo | ||
deallocate(Timelist) | ||
|
||
! swap around ref numbers for different data set | ||
ref_weights_leap(1) = day_before_leap_day | ||
ref_weights_leap(2) = day_before_leap_day ! feb 28th | ||
ref_weights_leap(3) = 0.0 | ||
ref_weights_leap(4) = day_before_leap_day_with_ly | ||
ref_weights_leap(5) = 0.0 | ||
ref_weights_leap(6) = 0.0 | ||
! Tests of modulo time and leap year inconsistency | ||
Time_beg = set_date(1978, 1, 1) | ||
Time_end = set_date(1981, 1, 1) | ||
|
@@ -210,6 +261,10 @@ program test_time_interp | |
call print_date(Time(ntest),' Time=') | ||
write(outunit,99) index1,index2,weight | ||
write(outunit,'()') | ||
if( .not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) & | ||
call mpp_error(FATAL, 'invalid results for indices with leap year correction') | ||
if(abs(weight - ref_weights_leap(ntest)) .gt. SMALL) & | ||
call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference") | ||
enddo | ||
deallocate(Timelist) | ||
|
||
|
@@ -297,4 +352,30 @@ subroutine diagram(nline,ntest,modulo_time) | |
|
||
end subroutine diagram | ||
|
||
end program test_time_interp | ||
!> helper function to check results | ||
!! true if invalid , false for valid | ||
logical function is_valid_indices(ind1, ind2, tList, tintv, res_weight, mtime) | ||
integer, intent(in) :: ind1, ind2 | ||
type(time_type), intent(in) :: tList(:), tintv | ||
real(TI_TEST_KIND_), intent(in) :: res_weight | ||
integer, intent(in) :: mtime | ||
integer :: i | ||
|
||
! modulo_time determines wrap around | ||
if( mtime .eq. NONE) then | ||
if (ind1 .eq. SIZE(tList)) then | ||
is_valid_indices = ind2 .eq. ind1 | ||
else | ||
is_valid_indices = ind2 .eq. ind1+1 | ||
endif | ||
else ! YEAR, default | ||
if (ind1 .eq. 12 ) then | ||
is_valid_indices = ind2 .eq. 1 | ||
else | ||
is_valid_indices = ind2 .eq. ind1+1 | ||
endif | ||
endif | ||
|
||
end function is_valid_indices | ||
|
||
end program test_time_interp |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ah it did that for you too! the do loops are repeated here
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
i think something wonky might be happening with our branches, might be better to simplify things and just use mixedmode as the starting point going forward