Skip to content

Commit

Permalink
bugfix in get_hit_times, l_correct_time
Browse files Browse the repository at this point in the history 
So far when specifying the sampling time for a dedicated output, like dt_movie,
the time array for the outputs was constructed assuming that the code would
run dt_old*n_time_steps, which could be quite wrong when the time step increases.

To prevent this, I rather fix the time array to one single entry in that case
and increment its value by dt_* whenever l_correct_time is hit.
  • Loading branch information
@tgastine
tgastine committed Jul 29, 2024
1 parent 9d930fa commit 8dce5a9
Show file tree
Hide file tree
Showing 5 changed files with 101 additions and 131 deletions.
2 changes: 1 addition & 1 deletion src/magic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,7 +392,7 @@ program magic
call initialize_courant(time, tscheme%dt(1), tag)

!--- Second pre-calculation:
call preCalcTimes(time, tEND, tscheme%dt(1), n_time_step, n_time_steps)
call preCalcTimes(time, n_time_step)

!--- Write info to STDOUT and log-file:
if ( rank == 0 ) then
Expand Down
22 changes: 11 additions & 11 deletions src/output_data.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,17 +26,17 @@ module output_data
real(cp), public :: t_TOmovie_start,t_TOmovie_stop,dt_TOmovie
real(cp), public :: t_pot_start,t_pot_stop,dt_pot
real(cp), public :: t_probe_start,t_probe_stop,dt_probe
integer, public :: n_graph_step,n_graphs,n_t_graph
integer, public :: n_rst_step,n_rsts,n_t_rst,n_stores
integer, public :: n_log_step,n_logs,n_t_log
integer, public :: n_spec_step,n_specs,n_t_spec
integer, public :: n_cmb_step,n_cmbs,n_t_cmb
integer, public :: n_r_field_step,n_r_fields,n_t_r_field
integer, public :: n_movie_step,n_movie_frames,n_t_movie
integer, public :: n_TO_step,n_TOs,n_t_TO
integer, public :: n_TOmovie_step,n_TOmovie_frames,n_t_TOmovie
integer, public :: n_pot_step,n_pots,n_t_pot
integer, public :: n_probe_step,n_probe_out,n_t_probe
integer, public :: n_graph_step,n_graphs
integer, public :: n_rst_step,n_rsts,n_stores
integer, public :: n_log_step,n_logs
integer, public :: n_spec_step,n_specs
integer, public :: n_cmb_step,n_cmbs
integer, public :: n_r_field_step,n_r_fields
integer, public :: n_movie_step,n_movie_frames
integer, public :: n_TO_step,n_TOs
integer, public :: n_TOmovie_step,n_TOmovie_frames
integer, public :: n_pot_step,n_pots
integer, public :: n_probe_step,n_probe_out
integer, public, parameter :: n_time_hits=30 ! Maximum number of specific times for I/O in input namelist
real(cp), allocatable, public :: t_graph(:)
real(cp), allocatable, public :: t_rst(:)
Expand Down
129 changes: 46 additions & 83 deletions src/preCalculations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ subroutine preCalc(tscheme)
!----- Rotational time scale:
tScale=one/ek ! or ekScaled ? (not defined yet...)
if ( rank==0 ) then
print*, 'Warning: rotational timescale, be sure to set dtmax large enough !'
write(output_unit,*) 'Warning: rotational timescale, be sure to set dtmax large enough !'
end if
end if
if ( n_lScale == 0 ) then
Expand Down Expand Up @@ -795,91 +795,76 @@ subroutine preCalc(tscheme)
end if
end subroutine preCalc
!-------------------------------------------------------------------------------
subroutine preCalcTimes(time,tEND,dt,n_time_step,n_time_steps)
subroutine preCalcTimes(time,n_time_step)
!
! Precalc. after time, time and dt has been read from startfile.
!

!-- Output variables
real(cp), intent(inout) :: time ! Current time
real(cp), intent(in) :: tEND ! Final requested time
real(cp), intent(in) :: dt ! Time step size
integer, intent(inout) :: n_time_step ! Index of time loop
integer, intent(in) :: n_time_steps ! Number of iterations

!-- Local variables:
real(cp) :: time_span

!----- Set time step:
if ( l_reset_t ) then
time=0.0_cp
n_time_step=0
end if

if ( tEND > 0.0_cp ) then
time_span = tEND-time
else
time_span = dt*n_time_steps
end if

tmagcon=tmagcon+time

!-- Get output times:
call get_hit_times(t_graph,n_t_graph,t_graph_start,t_graph_stop,dt_graph, &
& n_graphs,n_graph_step,'graph',time,time_span,tScale)
call get_hit_times(t_graph,t_graph_start,t_graph_stop,dt_graph, &
& n_graphs,n_graph_step,'graph',time,tScale)

call get_hit_times(t_pot,n_t_pot,t_pot_start,t_pot_stop,dt_pot, &
& n_pots,n_pot_step,'pot',time,time_span,tScale)
call get_hit_times(t_pot,t_pot_start,t_pot_stop,dt_pot, &
& n_pots,n_pot_step,'pot',time,tScale)

call get_hit_times(t_rst,n_t_rst,t_rst_start,t_rst_stop,dt_rst, &
& n_rsts,n_rst_step,'rst',time,time_span,tScale)
call get_hit_times(t_rst,t_rst_start,t_rst_stop,dt_rst, &
& n_rsts,n_rst_step,'rst',time,tScale)

call get_hit_times(t_log,n_t_log,t_log_start,t_log_stop,dt_log, &
& n_logs,n_log_step,'log',time,time_span,tScale)
call get_hit_times(t_log,t_log_start,t_log_stop,dt_log, &
& n_logs,n_log_step,'log',time,tScale)

call get_hit_times(t_spec,n_t_spec,t_spec_start,t_spec_stop,dt_spec, &
& n_specs,n_spec_step,'spec',time,time_span,tScale)
call get_hit_times(t_spec,t_spec_start,t_spec_stop,dt_spec, &
& n_specs,n_spec_step,'spec',time,tScale)

if ( l_probe ) then
call get_hit_times(t_probe,n_t_probe,t_probe_start,t_probe_stop,dt_probe, &
& n_probe_out,n_probe_step,'probe',time,time_span,tScale)
call get_hit_times(t_probe,t_probe_start,t_probe_stop,dt_probe, &
& n_probe_out,n_probe_step,'probe',time,tScale)
end if

if ( l_cmb_field ) then
l_cmb_field=.false.
call get_hit_times(t_cmb,n_t_cmb,t_cmb_start,t_cmb_stop,dt_cmb,n_cmbs,&
& n_cmb_step,'cmb',time,time_span,tScale)
if ( n_cmbs > 0 .or. n_cmb_step > 0 .or. n_t_cmb > 0 ) l_cmb_field= .true.
call get_hit_times(t_cmb,t_cmb_start,t_cmb_stop,dt_cmb,n_cmbs,&
& n_cmb_step,'cmb',time,tScale)
end if
l_dt_cmb_field=l_dt_cmb_field .and. l_cmb_field

if ( l_r_field ) then
call get_hit_times(t_r_field,n_t_r_field,t_r_field_start,t_r_field_stop,&
call get_hit_times(t_r_field,t_r_field_start,t_r_field_stop, &
& dt_r_field,n_r_fields,n_r_field_step,'r_field',time, &
& time_span,tScale)
& tScale)
end iF

if ( l_movie ) then
call get_hit_times(t_movie,n_t_movie,t_movie_start,t_movie_stop,dt_movie,&
& n_movie_frames,n_movie_step,'movie',time,time_span,tScale)
call get_hit_times(t_movie,t_movie_start,t_movie_stop,dt_movie,&
& n_movie_frames,n_movie_step,'movie',time,tScale)
end if

if ( l_TO ) then
if ( n_TOs == 0 .and. n_t_TO == 0 ) n_TO_step=max(3,n_TO_step)
call get_hit_times(t_TO,n_t_TO,t_TO_start,t_TO_stop,dt_TO,n_TOs, &
& n_TO_step,'TO',time,time_span,tScale)
call get_hit_times(t_TO,t_TO_start,t_TO_stop,dt_TO,n_TOs, &
& n_TO_step,'TO',time,tScale)
end if

if ( l_TOmovie ) then
call get_hit_times(t_TOmovie,n_t_TOmovie,t_TOmovie_start,t_TOmovie_stop,&
call get_hit_times(t_TOmovie,t_TOmovie_start,t_TOmovie_stop, &
& dt_TOmovie,n_TOmovie_frames,n_TOmovie_step,'TOmovie',&
& time,time_span,tScale)
& time,tScale)
end if

end subroutine preCalcTimes
!-------------------------------------------------------------------------------
subroutine get_hit_times(t,n_t,t_start,t_stop,dt,n_tot,n_step,string,time,&
& time_span,tScale)
subroutine get_hit_times(t,t_start,t_stop,dt,n_tot,n_step,string,time,&
& tScale)
!
! This subroutine checks whether any specific times t(*) are given
! on input. If so, it returns their number n_r and sets l_t to true.
Expand All @@ -890,7 +875,6 @@ subroutine get_hit_times(t,n_t,t_start,t_stop,dt,n_tot,n_step,string,time,&
!-- Input variables:
real(cp), intent(in) :: time ! Time of start file
real(cp), intent(in) :: tScale ! Scale unit for time
real(cp), intent(in) :: time_span ! Time interval
character(len=*), intent(in) :: string
integer, intent(inout) :: n_tot ! Number of output (times) if no times defined
integer, intent(inout) :: n_step ! Ouput step in no. of time steps
Expand All @@ -899,19 +883,25 @@ subroutine get_hit_times(t,n_t,t_start,t_stop,dt,n_tot,n_step,string,time,&
real(cp), intent(inout) :: t_stop ! Stop time for output
real(cp), intent(inout) :: dt ! Time step for output

!-- Output variables
integer, intent(out) :: n_t ! No. of output times

!-- Local variables:
logical :: l_t
real(cp) :: tmp(size(t))
integer :: n, n_t_max
integer :: n, n_t_max, n_t

n_t_max = size(t)
t_start=t_start/tScale
t_stop =t_stop/tScale
dt =dt/tScale

if ( n_tot /= 0 .and. n_step /= 0 ) then
write(output_unit,*)
write(output_unit,*) '! You have to either provide the total'
write(output_unit,*) '! number or the step for output:'
write(output_unit,'(A,2(A,I10))') string, "n_tot = ",n_tot,", n_step = ",n_step
write(output_unit,*) '! I set the step width to zero!'
n_step=0
end if

!-- Check whether any time is given explicitly:
l_t=.false.
n_t=0
Expand All @@ -935,58 +925,31 @@ subroutine get_hit_times(t,n_t,t_start,t_stop,dt,n_tot,n_step,string,time,&
if ( t_start < time ) t_start=time
if ( .not. l_t .and. ( dt > 0.0_cp .or. ( n_tot > 0 .and. t_stop > t_start ) ) ) then
if ( n_tot > 0 .and. dt > 0.0_cp ) then
n_t =n_tot
n_tot=0
else if ( dt > 0.0_cp ) then
if ( t_stop > t_start ) then
n_t=int((t_stop-t_start)/dt)+1
else ! In case t_start, t_stop are not provided guess a max number
n_t=int(time_span/dt)+1
end if
n_t=n_tot
n_tot=0 ! This is to ensure that time array is used later in l_correct_step
else if ( dt > 0.0_cp ) then ! In that case t(:) is an array with one element
n_t=1
else if ( n_tot > 0 ) then
n_t=n_tot
n_tot=0
dt=(t_stop-t_start)/real(n_t-1,kind=cp)
n_tot=0 ! This is to ensure that time array is used later in l_correct_step
end if
l_t=.true.

deallocate(t) ! Deallocate to get the proper size
if ( t_start == time ) then
n_t=n_t-1
allocate(t(n_t))
t(1)=t_start+dt
else
allocate(t(n_t))
t(1)=t_start
end if

allocate(t(n_t))
t(1)=t_start
do n=2,n_t
t(n)=t(n-1)+dt
end do

end if

if ( n_tot /= 0 .and. n_step /= 0 ) then
write(output_unit,*)
write(output_unit,*) '! You have to either provide the total'
write(output_unit,*) '! number or the step for output:'
write(output_unit,'(A,2(A,I10))') string, "n_tot = ",n_tot,", n_step = ",n_step
write(output_unit,*) '! I set the step width to zero!'
n_step=0
end if

if ( l_t ) then
t_start=t(1)
t_stop =t(n_t)
if ( n_t >= 2 ) then
dt=t(2)-t(1)
else
dt=0.0_cp
end if
else ! If this is still filled with -1, reduce the size to one element
!-- In case only n_step or n_tot is specified such that l_t=.false.
if ( .not. l_t ) then
deallocate(t)
allocate(t(1))
t(:)=-one
t(1)=-one
end if

end subroutine get_hit_times
Expand Down
Loading

0 comments on commit 8dce5a9

Please sign in to comment.