From 8dce5a95023b2244638c65145dfc8a51097bb7ea Mon Sep 17 00:00:00 2001 From: Thomas Gastine Date: Mon, 29 Jul 2024 17:00:01 +0200 Subject: [PATCH] bugfix in get_hit_times, l_correct_time 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. --- src/magic.f90 | 2 +- src/output_data.f90 | 22 +++---- src/preCalculations.f90 | 129 ++++++++++++++-------------------------- src/step_time.f90 | 62 +++++++++---------- src/useful.f90 | 17 ++++-- 5 files changed, 101 insertions(+), 131 deletions(-) diff --git a/src/magic.f90 b/src/magic.f90 index a72098ac..22984c34 100644 --- a/src/magic.f90 +++ b/src/magic.f90 @@ -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 diff --git a/src/output_data.f90 b/src/output_data.f90 index 0948e033..9acf4117 100644 --- a/src/output_data.f90 +++ b/src/output_data.f90 @@ -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(:) diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index 1d05a0f0..74391df0 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -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 @@ -795,20 +795,14 @@ 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 @@ -816,70 +810,61 @@ subroutine preCalcTimes(time,tEND,dt,n_time_step,n_time_steps) 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. @@ -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 @@ -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 @@ -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 diff --git a/src/step_time.f90 b/src/step_time.f90 index 6c9d004d..5625d17d 100644 --- a/src/step_time.f90 +++ b/src/step_time.f90 @@ -39,17 +39,17 @@ module step_time_mod & assemble_stage_Rdist use signals_mod, only: initialize_signals, check_signals use graphOut_mod, only: open_graph_file, close_graph_file - use output_data, only: tag, n_graph_step, n_graphs, n_t_graph, t_graph, & - & n_spec_step, n_specs, n_t_spec, t_spec, & - & n_movie_step, n_movie_frames, n_t_movie, t_movie,& - & n_TOmovie_step, n_TOmovie_frames, n_t_TOmovie, & - & t_TOmovie, n_pot_step, n_pots, n_t_pot, t_pot, & - & n_rst_step, n_rsts, n_t_rst, t_rst, n_stores, & - & n_log_step, n_logs, n_t_log, t_log, n_cmb_step, & - & n_cmbs, n_t_cmb, t_cmb, n_r_field_step, & - & n_r_fields, n_t_r_field, t_r_field, n_TO_step, & - & n_TOs, n_t_TO, t_TO, n_probe_step, n_probe_out, & - & n_t_probe, t_probe, log_file, n_log_file, & + use output_data, only: tag, n_graph_step, n_graphs, dt_graph, t_graph, & + & n_spec_step, n_specs, dt_spec, t_spec, & + & n_movie_step, n_movie_frames, dt_movie, t_movie,& + & n_TOmovie_step, n_TOmovie_frames, dt_TOmovie, & + & t_TOmovie, n_pot_step, n_pots, dt_pot, t_pot, & + & n_rst_step, n_rsts, dt_rst, t_rst, n_stores, & + & n_log_step, n_logs, dt_log, t_log, n_cmb_step, & + & n_cmbs, dt_cmb, t_cmb, n_r_field_step, & + & n_r_fields, dt_r_field, t_r_field, n_TO_step, & + & n_TOs, dt_TO, t_TO, n_probe_step, n_probe_out, & + & dt_probe, t_probe, log_file, n_log_file, & & n_time_hits use updateB_mod, only: get_mag_rhs_imp, get_mag_ic_rhs_imp, b_ghost, aj_ghost, & & get_mag_rhs_imp_ghost, fill_ghosts_B @@ -299,17 +299,17 @@ subroutine step_time(time, tscheme, n_time_step, run_time_start) if ( time >= tEND .and. tEND /= 0.0_cp ) l_stop_time=.true. !-- Checking logic for output: - l_graph= l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_graph_step,n_graphs,n_t_graph,t_graph) .or. & + l_graph= l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & + & n_graph_step,n_graphs,dt_graph,t_graph) .or. & & n_graph_signal == 1 n_graph_signal=0 ! reset interrupt signal ! l_spectrum= & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_spec_step,n_specs,n_t_spec,t_spec) .or. & + & n_spec_step,n_specs,dt_spec,t_spec) .or. & & n_spec_signal == 1 - l_frame= l_movie .and. ( & - & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_movie_step,n_movie_frames,n_t_movie,t_movie) .or. & + l_frame= l_movie .and. ( & + & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & + & n_movie_step,n_movie_frames,dt_movie,t_movie) .or. & & n_time_steps_go == 1 ) if ( l_mag .or. l_mag_LF ) then l_dtB=( l_frame .and. l_dtBmovie ) .or. ( l_log .and. l_DTrMagSpec ) @@ -317,54 +317,54 @@ subroutine step_time(time, tscheme, n_time_step, run_time_start) lTOframe=l_TOmovie .and. & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_TOmovie_step,n_TOmovie_frames,n_t_TOmovie,t_TOmovie) + & n_TOmovie_step,n_TOmovie_frames,dt_TOmovie,t_TOmovie) l_probe_out=l_probe .and. & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_probe_step,n_probe_out,n_t_probe,t_probe) + & n_probe_step,n_probe_out,dt_probe,t_probe) !-- Potential files l_pot= l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_pot_step,n_pots,n_t_pot,t_pot) .or. & + & n_pot_step,n_pots,dt_pot,t_pot) .or. & & n_pot_signal == 1 n_pot_signal=0 ! reset interrupt signal ! l_new_rst_file= & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_rst_step,n_rsts,n_t_rst,t_rst) .or. & + & n_rst_step,n_rsts,dt_rst,t_rst) .or. & & n_rst_signal == 1 n_rst_signal=0 l_store= l_new_rst_file .or. & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & 0,n_stores,0,t_rst) + & 0,n_stores,0.0_cp,t_rst) l_log= l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_log_step,n_logs,n_t_log,t_log) + & n_log_step,n_logs,dt_log,t_log) l_cmb= l_cmb_field .and. & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_cmb_step,n_cmbs,n_t_cmb,t_cmb) + & n_cmb_step,n_cmbs,dt_cmb,t_cmb) l_r= l_r_field .and. & & l_correct_step(n_time_step-1,time,timeLast,n_time_steps, & - & n_r_field_step,n_r_fields,n_t_r_field, & + & n_r_field_step,n_r_fields,dt_r_field, & & t_r_field) l_logNext=.false. if ( n_time_step+1 <= n_time_steps+1 ) & & l_logNext= & & l_correct_step(n_time_step,time+tscheme%dt(1),timeLast, & - & n_time_steps,n_log_step,n_logs,n_t_log,t_log) + & n_time_steps,n_log_step,n_logs,dt_log,t_log) lTOCalc= n_time_step > 2 .and. l_TO .and. & & l_correct_step(n_time_step-1,time,timeLast, & - & n_time_steps,n_TO_step,n_TOs,n_t_TO,t_TO) + & n_time_steps,n_TO_step,n_TOs,dt_TO,t_TO) lTOnext =.false. lTOframeNext=.false. if ( n_time_step+1 <= n_time_steps+1 ) then lTONext= l_TO .and. & & l_correct_step(n_time_step,time+tscheme%dt(1),& - & timeLast,n_time_steps,n_TO_step,n_TOs,n_t_TO,t_TO) + & timeLast,n_time_steps,n_TO_step,n_TOs,dt_TO,t_TO) lTOframeNext= l_TOmovie .and. & & l_correct_step(n_time_step,time+tscheme%dt(1), & & timeLast,n_time_steps,n_TOmovie_step, & - & n_TOmovie_frames,n_t_TOmovie,t_TOmovie) + & n_TOmovie_frames,dt_TOmovie,t_TOmovie) end if lTONext =lTOnext.or.lTOframeNext lTONext2 =.false. @@ -373,11 +373,11 @@ subroutine step_time(time, tscheme, n_time_step, run_time_start) lTONext2= l_TO .and. & & l_correct_step(n_time_step+1,time+2*tscheme%dt(1), & & timeLast,n_time_steps,n_TO_step, & - & n_TOs,n_t_TO,t_TO) + & n_TOs,dt_TO,t_TO) lTOframeNext2= l_TOmovie .and. & & l_correct_step(n_time_step+1,time+2*tscheme%dt(1), & & timeLast,n_time_steps,n_TOmovie_step, & - & n_TOmovie_frames,n_t_TOmovie,t_TOmovie) + & n_TOmovie_frames,dt_TOmovie,t_TOmovie) end if lTONext2=lTOnext2.or.lTOframeNext2 diff --git a/src/useful.f90 b/src/useful.f90 index 1965660e..f7114390 100644 --- a/src/useful.f90 +++ b/src/useful.f90 @@ -20,7 +20,7 @@ module useful contains logical function l_correct_step(n,t,t_last,n_max,n_step,n_intervals, & - & n_ts,times) + & dt_output,times) ! ! Suppose we have a (loop) maximum of n_max steps! ! If n_intervals times in these steps a certain action should be carried out @@ -41,8 +41,8 @@ logical function l_correct_step(n,t,t_last,n_max,n_step,n_intervals, & integer, intent(in) :: n_max ! max number of steps integer, intent(in) :: n_step ! action interval integer, intent(in) :: n_intervals ! number of actions - integer, intent(in) :: n_ts ! number of times t - real(cp), intent(in) :: times(:) ! times where l_correct_step == true + real(cp), intent(in) :: dt_output ! Output frequency + real(cp), intent(inout) :: times(:) ! times where l_correct_step == true !-- Local variables: integer :: n_offset ! offset with no action @@ -66,8 +66,15 @@ logical function l_correct_step(n,t,t_last,n_max,n_step,n_intervals, & if ( n == n_max .or. mod(n,n_step) == 0 ) l_correct_step=.true. end if - if ( n_ts >= 1 ) then - do n_t=1,n_ts + if ( size(times) == 1 .and. times(1) > 0.0_cp ) then + !-- Time array has one single entry for the next output + if ( times(1) < t .and. times(1) >= t_last ) then + l_correct_step=.true. + times(1)=times(1)+dt_output ! In this case, increment the target time for next output + end if + else if ( size(times) > 1 ) then + !-- Time array contains multiple entries + do n_t=1,size(times) if ( times(n_t) < t .and. times(n_t) >= t_last ) then l_correct_step=.true. exit