Skip to content
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.

Sign up for GitHub

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

Evaluate time dependent tendencies at correct RK4 stage time value #55

Closed
wants to merge 10 commits into from
Closed

Evaluate time dependent tendencies at correct RK4 stage time value #55

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
b9062f1
Add RK4 stage time evaluation factors
sbrus89 Jun 9, 2023
1d28e13
Add RK4 stage time evaluation to manufactured solution
sbrus89 Jun 9, 2023
4039c01
Add RK stage time evaluation to tidal potential
sbrus89 Jun 9, 2023
710a369
Add RK4 stage time evaluation to time varying forcing
sbrus89 Jun 9, 2023
6456826
Fix unassociated forcingHead when time varying forcing is off
sbrus89 Jun 13, 2023
e663090
Fix ordering of RK4 time varying forcing call
sbrus89 Jun 22, 2023
3b91be0
Apply fix in mpas_set_timeInterval for negative fractional timesteps
sbrus89 Jul 3, 2023
9d21c15
Add _RKIND to RK4 weights
sbrus89 Oct 14, 2024
7773fe4
Remove unneeded decimal timestep code
sbrus89 Oct 14, 2024
c62698d
Use enum for time integrator choice
sbrus89 Oct 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 45 additions & 6 deletions components/mpas-framework/src/framework/mpas_forcing.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@ module mpas_forcing
mpas_forcing_init_field_data, &
mpas_forcing_get_forcing, &
mpas_forcing_get_forcing_time, &
mpas_forcing_write_restart_times
mpas_forcing_write_restart_times, &
mpas_advance_forcing_clock

contains

Expand Down Expand Up @@ -1193,7 +1194,7 @@ subroutine mpas_forcing_get_forcing(&!{{{
if (trim(forcingGroup % forcingGroupName) == trim(forcingGroupName)) then

! advance the forcing time
call advance_forcing_clock(forcingGroup, dt)
call mpas_advance_forcing_clock(forcingGroup, dt)

! cycle the forcing clock
if (forcingGroup % forcingCycleUse) then
Expand Down Expand Up @@ -1230,7 +1231,7 @@ end subroutine mpas_forcing_get_forcing!}}}

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! advance_forcing_clock
! mpas_advance_forcing_clock
!
!> \brief set the forcing clock
!> \author Adrian K. Turner, LANL
Expand All @@ -1241,7 +1242,7 @@ end subroutine mpas_forcing_get_forcing!}}}
!
!-----------------------------------------------------------------------

subroutine advance_forcing_clock(&!{{{
subroutine mpas_advance_forcing_clock(&!{{{
forcingGroup, &
dt)

Expand All @@ -1254,11 +1255,47 @@ subroutine advance_forcing_clock(&!{{{
type(MPAS_TimeInterval_type) :: &
timeStep ! time step interval

!! assuming it is not possible to give dts of months or years
!integer :: DD, H, M, S, S_n, S_d, foundNumAndDen, powerOfTen

!real(kind=RKIND) :: factor

!DD = dt / 86400_RKIND
!H = dt / 3600_RKIND
!M = dt / 60_RKIND
!S = dt
!S_n = 0
!S_d = 0
!if (abs(real(S) - dt) > 1.e-10) then !the time step has decimals
! S = 0
! foundNumAndDen = 0
! powerOfTen = 1
! factor = 10_RKIND
! S_n = abs(dt) * factor
! S_d = factor
! do while (foundNumAndDen == 0 .and. powerOfTen < 11)
! if (abs(real(S_n)/real(S_d) - abs(dt)) < 1.e-10) then
! foundNumAndDen = 1
! else
! powerOfTen = powerOfTen + 1
! factor = 10_RKIND ** powerOfTen
! S_n = abs(dt) * factor
! S_d = factor
! end if
! end do
! if (dt < 0.0_RKIND) then
! S_n = - S_n
! end if
!end if

sbrus89 marked this conversation as resolved.
Show resolved Hide resolved
! increment clock with timestep
! call mpas_set_timeInterval(timeStep, DD=DD, H=H, M=M, S=S, S_n=S_n, S_d=S_d)
call mpas_set_timeInterval(timeStep, dt=dt)
!call mpas_log_write('jl time step $i $i $i', intArgs=(/S, S_n, S_d /))
!call mpas_log_write('time step $i $i $i', intArgs=(/int(timeStep%ti%basetime%S), int(timeStep%ti%basetime%Sn), int(timeStep%ti%basetime%Sd) /))
call mpas_advance_clock(forcingGroup % forcingClock, timeStep)

end subroutine advance_forcing_clock!}}}
end subroutine mpas_advance_forcing_clock!}}}

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
Expand Down Expand Up @@ -1597,8 +1634,10 @@ subroutine get_interpolants_linear(interpolants, forcingStream, currentTime)!{{{
call mpas_get_timeInterval(diff1, forcingStream % forcingTimes(1), dt=diffr1)
call mpas_get_timeInterval(diff2, currentTime, dt=diffr2)

!call mpas_log_write('diffr2 $r, diffr, $r', realArgs=(/ diffr2, diffr /))
interpolants(1) = diffr2 / diffr
interpolants(2) = 1.0_RKIND - interpolants(1) !diffr1 / diffr


end subroutine get_interpolants_linear!}}}

Expand All @@ -1612,7 +1651,7 @@ end subroutine get_interpolants_linear!}}}
!> \details
!> Given the current time and forcing times calculate the correct
!> interpolants with piecewise constant interpolation
!

!-----------------------------------------------------------------------

subroutine get_interpolants_constant(interpolants, forcingStream, currentTime)!{{{
Expand Down
8 changes: 8 additions & 0 deletions components/mpas-framework/src/framework/mpas_timekeeping.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1427,6 +1427,8 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8,
timeString_ = timeString
end if

!call mpas_log_write('timeString_ '//trim(timeString_))

numerator = 0
denominator = 1

Expand Down Expand Up @@ -1519,6 +1521,12 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8,
return
end if

if (sec < 0) then
numerator = -numerator
end if

!call mpas_log_write('sec $i, num $i, denom $i', intArgs=(/int(sec),numerator,denominator/))

call ESMF_TimeIntervalSet(interval % ti, YY=year, MM=month, D=day, H=hour, M=min, S_i8=sec, Sn=numerator, Sd=denominator, rc=ierr)

else
Expand Down
4 changes: 4 additions & 0 deletions components/mpas-ocean/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4130,6 +4130,10 @@
description="Multiplication factors to smooth ssh at coastlines for SAL caculation"
packages="tidalPotentialForcingPKG"
/>
<var name="forcingTimeIncrement" type="real" dimensions="" units="s"
description="Time increment to compute forcing terms at intermediate time intervals"
packages="forwardMode"
/>
</var_struct>
<var_struct name="timeVaryingForcing" time_levs="1">
<var name="windSpeedU" type="real" dimensions="nCells Time" units="m s^-1"
Expand Down
21 changes: 19 additions & 2 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ module ocn_forward_mode
use mpas_stream_manager
use mpas_timekeeping
use mpas_dmpar
use mpas_forcing
use mpas_timer
use mpas_log
use mpas_decomp
Expand Down Expand Up @@ -94,6 +95,7 @@ module ocn_forward_mode

use ocn_forcing
use ocn_time_varying_forcing
use ocn_framework_forcing

use ocn_constants
use ocn_config
Expand Down Expand Up @@ -664,7 +666,12 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{

! initialize time-varying forcing
call ocn_time_varying_forcing_init(domain)
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)

! if not using RK4, calculate time varying forcing terms once per
! time-step as opposed at each RK substage as implemented in RK4
if (timeIntegratorChoice /= 4) then
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if now that we have the capability to compute the forcing correctly, we should allow the possibility of not doing it. For instance, instead of having a condition on the type of integrator here, we could have something like config_update_forcing_at_intermediate_timesteps that is also used in the RK4 timestepping code and sets all the forcingTimeIncrementRK4 to zero if false. In this way, for coupled runs, one can decide to have a lighter coupling, but if one wants to do convergence tests and use RK4 as a reference solution to validate other time stepping schemes it is still possible turning that flag on.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (timeIntegratorChoice /= 4) then
if (timeIntegratorChoice /= timeIntRK4) then

call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
endif

! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
Expand Down Expand Up @@ -834,7 +841,17 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
endif

! read in next time level data required for time-varying forcing
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
if (timeIntegratorChoice /= 4) then
sbrus89 marked this conversation as resolved.
Show resolved Hide resolved
! if not using RK4, calculate time varying forcing terms once per
! time-step as opposed at each RK substage as implemented in RK4
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
else
if (config_use_time_varying_atmospheric_forcing .or. &
config_use_time_varying_land_ice_forcing) then
! increment forcing clock to next time-step
call mpas_advance_forcing_clock(forcingGroupHead, dt)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, maybe instead of setting the forcingTimeIncrementRK4 to zero in terms of having a light coupling we should not advance the forcing clock here. Please check if this statement is true. Thanks!

endif
endif

! Validate that the state is OK to run with for the next timestep.
call ocn_validate_state(domain, timeLevel=1)
Expand Down
6 changes: 3 additions & 3 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,9 @@ module ocn_time_integration
!
!--------------------------------------------------------------------

! Enum for selecting different time integrators
integer, public :: timeIntegratorChoice

!--------------------------------------------------------------------
!
! Public member functions
Expand All @@ -63,9 +66,6 @@ module ocn_time_integration
!
!--------------------------------------------------------------------

! Enum for selecting different time integrators
integer :: timeIntegratorChoice

integer, parameter :: &
timeIntUnknown = 0, &! unknown or undefined
timeIntSplitExplicit = 1, &! split-explicit
sbrus89 marked this conversation as resolved.
Show resolved Hide resolved
Expand Down
27 changes: 25 additions & 2 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ module ocn_time_integration_rk4
use ocn_effective_density_in_land_ice
use ocn_surface_land_ice_fluxes
use ocn_transport_tests
use ocn_time_varying_forcing

use ocn_subgrid

Expand Down Expand Up @@ -118,7 +119,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{

type (mpas_pool_type), pointer :: nextProvisPool, prevProvisPool

real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights, forcingTimeIncrementRK4

real (kind=RKIND) :: coef
real (kind=RKIND), dimension(:,:), pointer :: &
Expand Down Expand Up @@ -186,6 +187,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
! Tidal boundary condition
logical, pointer :: config_use_tidal_forcing
character (len=StrKIND), pointer :: config_tidal_forcing_type
real (kind=RKIND), pointer :: forcingTimeIncrement
real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue
real (kind=RKIND), dimension(:,:), pointer :: restingThickness
real (kind=RKIND) :: totalDepth
Expand Down Expand Up @@ -233,6 +235,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)

call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool)

Expand Down Expand Up @@ -260,6 +263,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
subgridSshCellTableRange)

call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)

forcingTimeIncrement = 0.0_RKIND

! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur
! assignments below are needed to maintain bit-for-bit results
Expand Down Expand Up @@ -360,6 +366,14 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
rk_substep_weights(3) = dt
rk_substep_weights(4) = dt ! a_4 only used for ALE step, otherwise it is skipped.

! these are time increments to evaluate the tidal forcing at the
! intermediate time-steps as required by RK4

forcingTimeIncrementRK4(1) = 0.0_RKIND
forcingTimeIncrementRK4(2) = dt/2.
forcingTimeIncrementRK4(3) = dt/2.
forcingTimeIncrementRK4(4) = dt

call mpas_timer_start("RK4-main loop")

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Expand Down Expand Up @@ -481,8 +495,17 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{

block => domain % blocklist
do while (associated(block))
call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
forcingTimeIncrement = forcingTimeIncrementRK4(rk_step)
block => block % next
end do

call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)

block => domain % blocklist
do while (associated(block))
call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err )
block => block % next
end do
Expand Down
7 changes: 7 additions & 0 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -380,6 +380,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{

real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew

real (kind=RKIND), pointer :: forcingTimeIncrement

! Remap variables
real (kind=RKIND), dimension(:,:), pointer :: &
layerThicknessLagNew
Expand Down Expand Up @@ -507,6 +509,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{

call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2)

call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
forcingTimeIncrement)

forcingTimeIncrement = 0.0_RKIND

allocate(bottomDepthEdge(nEdgesAll+1))

if (config_transport_tests_flow_id > 0) then
Expand Down
7 changes: 7 additions & 0 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,6 +280,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew

real (kind=RKIND), pointer :: forcingTimeIncrement

! Remap variables
real (kind=RKIND), dimension(:,:), pointer :: &
layerThicknessLagNew
Expand Down Expand Up @@ -398,6 +400,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

allocate(baroclinicThickness(nEdgesAll+1))

call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
forcingTimeIncrement)

forcingTimeIncrement = 0.0_RKIND

if (config_transport_tests_flow_id > 0) then
! This is a transport test. Write advection velocity from prescribed
! flow field.
Expand Down
26 changes: 22 additions & 4 deletions components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,13 @@ module ocn_manufactured_solution
!
!-----------------------------------------------------------------------

subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{
subroutine ocn_manufactured_solution_tend_thick(forcingPool, tend, err)!{{{

!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------

type (mpas_pool_type), intent(in) :: forcingPool

!-----------------------------------------------------------------
! input/output variables
Expand All @@ -97,14 +103,17 @@ subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{

integer :: iCell, kmin, kmax, k
real (kind=RKIND) :: phase, time
real (kind=RKIND), pointer :: forcingTimeIncrement

! End preamble
!-------------
! Begin code

if (.not. config_use_manufactured_solution) return

time = daysSinceStartOfSim*86400.0_RKIND
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)

time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement

do iCell = 1,nCellsOwned

Expand Down Expand Up @@ -137,7 +146,13 @@ end subroutine ocn_manufactured_solution_tend_thick!}}}
!> This routine computes the velocity tendency for the manufactured solution
!
!-----------------------------------------------------------------------
subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{
subroutine ocn_manufactured_solution_tend_vel(forcingPool, tend, err)!{{{

!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------

type (mpas_pool_type), intent(in) :: forcingPool

!-----------------------------------------------------------------
! input/output variables
Expand All @@ -158,14 +173,17 @@ subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{

integer :: iEdge, kmin, kmax, k
real (kind=RKIND) :: phase, u, v, time
real (kind=RKIND), pointer :: forcingTimeIncrement

! End preamble
!-----------------------------------------------------------------
! Begin code

if (.not. config_use_manufactured_solution) return

time = daysSinceStartOfSim*86400.0_RKIND
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)

time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement

do iEdge = 1, nEdgesOwned

Expand Down
Loading