Skip to content

Commit b77743a

Browse files
jonbobjonbob
committed
Merge branch 'sbrus89/ocn/rk4-time-fix' into next (PR #6695)
Evaluate time-dependent tendencies at correct time for each RK4 stage This PR corrects a longstanding issue with RK4 and tendency terms which have an explicit time dependence. Big thanks to @gcapodag and @jeremy-lilly for providing the fix in their LTS development branch. Here, I have migrated these changes from their MPAS-Dev branch to E3SM and added the modifications to correct the convergence issues for the manufactured solution test case [BFB] - mpaso standalone
2 parents c84b9d7 + c62698d commit b77743a

12 files changed

+194
-46
lines changed

‎components/mpas-framework/src/framework/mpas_forcing.F‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,8 @@ module mpas_forcing
3535
mpas_forcing_init_field_data, &
3636
mpas_forcing_get_forcing, &
3737
mpas_forcing_get_forcing_time, &
38-
mpas_forcing_write_restart_times
38+
mpas_forcing_write_restart_times, &
39+
mpas_advance_forcing_clock
3940

4041
contains
4142

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

11951196
! advance the forcing time
1196-
call advance_forcing_clock(forcingGroup, dt)
1197+
call mpas_advance_forcing_clock(forcingGroup, dt)
11971198

11981199
! cycle the forcing clock
11991200
if (forcingGroup % forcingCycleUse) then
@@ -1230,7 +1231,7 @@ end subroutine mpas_forcing_get_forcing!}}}
12301231

12311232
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
12321233
!
1233-
! advance_forcing_clock
1234+
! mpas_advance_forcing_clock
12341235
!
12351236
!> \brief set the forcing clock
12361237
!> \author Adrian K. Turner, LANL
@@ -1241,7 +1242,7 @@ end subroutine mpas_forcing_get_forcing!}}}
12411242
!
12421243
!-----------------------------------------------------------------------
12431244

1244-
subroutine advance_forcing_clock(&!{{{
1245+
subroutine mpas_advance_forcing_clock(&!{{{
12451246
forcingGroup, &
12461247
dt)
12471248

@@ -1256,9 +1257,10 @@ subroutine advance_forcing_clock(&!{{{
12561257

12571258
! increment clock with timestep
12581259
call mpas_set_timeInterval(timeStep, dt=dt)
1260+
!call mpas_log_write('time step from mpas_set_timeInterval: $i $i $i', intArgs=(/int(timeStep%ti%basetime%S), int(timeStep%ti%basetime%Sn), int(timeStep%ti%basetime%Sd) /))
12591261
call mpas_advance_clock(forcingGroup % forcingClock, timeStep)
12601262

1261-
end subroutine advance_forcing_clock!}}}
1263+
end subroutine mpas_advance_forcing_clock!}}}
12621264

12631265
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
12641266
!
@@ -1597,8 +1599,10 @@ subroutine get_interpolants_linear(interpolants, forcingStream, currentTime)!{{{
15971599
call mpas_get_timeInterval(diff1, forcingStream % forcingTimes(1), dt=diffr1)
15981600
call mpas_get_timeInterval(diff2, currentTime, dt=diffr2)
15991601

1602+
!call mpas_log_write('diffr2 $r, diffr, $r', realArgs=(/ diffr2, diffr /))
16001603
interpolants(1) = diffr2 / diffr
16011604
interpolants(2) = 1.0_RKIND - interpolants(1) !diffr1 / diffr
1605+
16021606

16031607
end subroutine get_interpolants_linear!}}}
16041608

@@ -1612,7 +1616,7 @@ end subroutine get_interpolants_linear!}}}
16121616
!> \details
16131617
!> Given the current time and forcing times calculate the correct
16141618
!> interpolants with piecewise constant interpolation
1615-
!
1619+
16161620
!-----------------------------------------------------------------------
16171621

16181622
subroutine get_interpolants_constant(interpolants, forcingStream, currentTime)!{{{

‎components/mpas-framework/src/framework/mpas_timekeeping.F‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1427,6 +1427,8 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8,
14271427
timeString_ = timeString
14281428
end if
14291429

1430+
!call mpas_log_write('timeString_ '//trim(timeString_))
1431+
14301432
numerator = 0
14311433
denominator = 1
14321434

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

1524+
if (sec < 0) then
1525+
numerator = -numerator
1526+
end if
1527+
1528+
!call mpas_log_write('sec $i, num $i, denom $i', intArgs=(/int(sec),numerator,denominator/))
1529+
15221530
call ESMF_TimeIntervalSet(interval % ti, YY=year, MM=month, D=day, H=hour, M=min, S_i8=sec, Sn=numerator, Sd=denominator, rc=ierr)
15231531

15241532
else

‎components/mpas-ocean/src/Registry.xml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4175,6 +4175,10 @@
41754175
description="Multiplication factors to smooth ssh at coastlines for SAL caculation"
41764176
packages="tidalPotentialForcingPKG"
41774177
/>
4178+
<var name="forcingTimeIncrement" type="real" dimensions="" units="s"
4179+
description="Time increment to compute forcing terms at intermediate time intervals"
4180+
packages="forwardMode"
4181+
/>
41784182
</var_struct>
41794183
<var_struct name="timeVaryingForcing" time_levs="1">
41804184
<var name="windSpeedU" type="real" dimensions="nCells Time" units="m s^-1"

‎components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F‎

Lines changed: 19 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ module ocn_forward_mode
2626
use mpas_stream_manager
2727
use mpas_timekeeping
2828
use mpas_dmpar
29+
use mpas_forcing
2930
use mpas_timer
3031
use mpas_log
3132
use mpas_decomp
@@ -94,6 +95,7 @@ module ocn_forward_mode
9495

9596
use ocn_forcing
9697
use ocn_time_varying_forcing
98+
use ocn_framework_forcing
9799

98100
use ocn_constants
99101
use ocn_config
@@ -664,7 +666,12 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
664666

665667
! initialize time-varying forcing
666668
call ocn_time_varying_forcing_init(domain)
667-
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
669+
670+
! if not using RK4, calculate time varying forcing terms once per
671+
! time-step as opposed at each RK substage as implemented in RK4
672+
if (timeIntegratorChoice /= timeIntRK4) then
673+
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
674+
endif
668675

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

836843
! read in next time level data required for time-varying forcing
837-
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
844+
if (timeIntegratorChoice /= timeIntRK4) then
845+
! if not using RK4, calculate time varying forcing terms once per
846+
! time-step as opposed at each RK substage as implemented in RK4
847+
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
848+
else
849+
if (config_use_time_varying_atmospheric_forcing .or. &
850+
config_use_time_varying_land_ice_forcing) then
851+
! increment forcing clock to next time-step
852+
call mpas_advance_forcing_clock(forcingGroupHead, dt)
853+
endif
854+
endif
838855

839856
! Validate that the state is OK to run with for the next timestep.
840857
call ocn_validate_state(domain, timeLevel=1)

‎components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F‎

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,19 @@ module ocn_time_integration
4848
!
4949
!--------------------------------------------------------------------
5050

51+
! Enum for selecting different time integrators
52+
integer, public :: timeIntegratorChoice
53+
54+
integer, public, parameter :: &
55+
timeIntUnknown = 0, &! unknown or undefined
56+
timeIntSplitExplicit = 1, &! split-explicit
57+
timeIntUnsplitExplicit = 2, &! unsplit-explicit
58+
timeIntSemiImplicit = 3, &! Semi-implicit
59+
timeIntRK4 = 4, &! 4th-order Runge-Kutta
60+
timeIntLTS = 5, &! local time-stepping
61+
timeIntFBLTS = 6, &! forward-backward lts
62+
timeIntSplitExplicitAB2 = 7 ! split-explicit AB2 baroclinic
63+
5164
!--------------------------------------------------------------------
5265
!
5366
! Public member functions
@@ -63,18 +76,6 @@ module ocn_time_integration
6376
!
6477
!--------------------------------------------------------------------
6578

66-
! Enum for selecting different time integrators
67-
integer :: timeIntegratorChoice
68-
69-
integer, parameter :: &
70-
timeIntUnknown = 0, &! unknown or undefined
71-
timeIntSplitExplicit = 1, &! split-explicit
72-
timeIntUnsplitExplicit = 2, &! unsplit-explicit
73-
timeIntSemiImplicit = 3, &! Semi-implicit
74-
timeIntRK4 = 4, &! 4th-order Runge-Kutta
75-
timeIntLTS = 5, &! local time-stepping
76-
timeIntFBLTS = 6, &! forward-backward lts
77-
timeIntSplitExplicitAB2 = 7 ! split-explicit AB2 baroclinic
7879

7980
!***********************************************************************
8081

‎components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F‎

Lines changed: 31 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ module ocn_time_integration_rk4
4343
use ocn_effective_density_in_land_ice
4444
use ocn_surface_land_ice_fluxes
4545
use ocn_transport_tests
46+
use ocn_time_varying_forcing
4647

4748
use ocn_subgrid
4849

@@ -118,7 +119,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
118119

119120
type (mpas_pool_type), pointer :: nextProvisPool, prevProvisPool
120121

121-
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
122+
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights, forcingTimeIncrementRK4
122123

123124
real (kind=RKIND) :: coef
124125
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -186,6 +187,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
186187
! Tidal boundary condition
187188
logical, pointer :: config_use_tidal_forcing
188189
character (len=StrKIND), pointer :: config_tidal_forcing_type
190+
real (kind=RKIND), pointer :: forcingTimeIncrement
189191
real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue
190192
real (kind=RKIND), dimension(:,:), pointer :: restingThickness
191193
real (kind=RKIND) :: totalDepth
@@ -233,6 +235,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
233235
call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
234236
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
235237
call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
238+
call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
236239

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

@@ -260,6 +263,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
260263
call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
261264
subgridSshCellTableRange)
262265

266+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
267+
268+
forcingTimeIncrement = 0.0_RKIND
263269

264270
! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur
265271
! assignments below are needed to maintain bit-for-bit results
@@ -343,10 +349,10 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
343349
! The coefficients of k_j are b_j = (1/6, 1/3, 1/3, 1/6) and are
344350
! initialized here as delta t * b_j:
345351

346-
rk_weights(1) = dt/6.
347-
rk_weights(2) = dt/3.
348-
rk_weights(3) = dt/3.
349-
rk_weights(4) = dt/6.
352+
rk_weights(1) = dt/6.0_RKIND
353+
rk_weights(2) = dt/3.0_RKIND
354+
rk_weights(3) = dt/3.0_RKIND
355+
rk_weights(4) = dt/6.0_RKIND
350356

351357
! The a_j coefficients of h in the computation of k_j are typically written (0, 1/2, 1/2, 1).
352358
! However, in the algorithm below we pre-compute the state for the tendency one iteration early.
@@ -355,11 +361,19 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
355361
! That is why the coefficients of h are one index early in the following, i.e.
356362
! a = (1/2, 1/2, 1)
357363

358-
rk_substep_weights(1) = dt/2.
359-
rk_substep_weights(2) = dt/2.
364+
rk_substep_weights(1) = dt/2.0_RKIND
365+
rk_substep_weights(2) = dt/2.0_RKIND
360366
rk_substep_weights(3) = dt
361367
rk_substep_weights(4) = dt ! a_4 only used for ALE step, otherwise it is skipped.
362368

369+
! these are time increments to evaluate the tidal forcing at the
370+
! intermediate time-steps as required by RK4
371+
372+
forcingTimeIncrementRK4(1) = 0.0_RKIND
373+
forcingTimeIncrementRK4(2) = dt/2.0_RKIND
374+
forcingTimeIncrementRK4(3) = dt/2.0_RKIND
375+
forcingTimeIncrementRK4(4) = dt
376+
363377
call mpas_timer_start("RK4-main loop")
364378

365379
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -481,8 +495,17 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
481495

482496
block => domain % blocklist
483497
do while (associated(block))
484-
call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
498+
call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
499+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
500+
forcingTimeIncrement = forcingTimeIncrementRK4(rk_step)
501+
block => block % next
502+
end do
503+
504+
call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
485505

506+
block => domain % blocklist
507+
do while (associated(block))
508+
call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
486509
call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err )
487510
block => block % next
488511
end do

‎components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -380,6 +380,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
380380

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

383+
real (kind=RKIND), pointer :: forcingTimeIncrement
384+
383385
! Remap variables
384386
real (kind=RKIND), dimension(:,:), pointer :: &
385387
layerThicknessLagNew
@@ -507,6 +509,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
507509

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

512+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
513+
forcingTimeIncrement)
514+
515+
forcingTimeIncrement = 0.0_RKIND
516+
510517
allocate(bottomDepthEdge(nEdgesAll+1))
511518

512519
if (config_transport_tests_flow_id > 0) then

‎components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -280,6 +280,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
280280

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

283+
real (kind=RKIND), pointer :: forcingTimeIncrement
284+
283285
! Remap variables
284286
real (kind=RKIND), dimension(:,:), pointer :: &
285287
layerThicknessLagNew
@@ -398,6 +400,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
398400

399401
allocate(baroclinicThickness(nEdgesAll+1))
400402

403+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
404+
forcingTimeIncrement)
405+
406+
forcingTimeIncrement = 0.0_RKIND
407+
401408
if (config_transport_tests_flow_id > 0) then
402409
! This is a transport test. Write advection velocity from prescribed
403410
! flow field.

‎components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F‎

Lines changed: 22 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,13 @@ module ocn_manufactured_solution
7676
!
7777
!-----------------------------------------------------------------------
7878

79-
subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{
79+
subroutine ocn_manufactured_solution_tend_thick(forcingPool, tend, err)!{{{
80+
81+
!-----------------------------------------------------------------
82+
! input variables
83+
!-----------------------------------------------------------------
84+
85+
type (mpas_pool_type), intent(in) :: forcingPool
8086

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

98104
integer :: iCell, kmin, kmax, k
99105
real (kind=RKIND) :: phase, time
106+
real (kind=RKIND), pointer :: forcingTimeIncrement
100107

101108
! End preamble
102109
!-------------
103110
! Begin code
104111

105112
if (.not. config_use_manufactured_solution) return
106113

107-
time = daysSinceStartOfSim*86400.0_RKIND
114+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
115+
116+
time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement
108117

109118
do iCell = 1,nCellsOwned
110119

@@ -137,7 +146,13 @@ end subroutine ocn_manufactured_solution_tend_thick!}}}
137146
!> This routine computes the velocity tendency for the manufactured solution
138147
!
139148
!-----------------------------------------------------------------------
140-
subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{
149+
subroutine ocn_manufactured_solution_tend_vel(forcingPool, tend, err)!{{{
150+
151+
!-----------------------------------------------------------------
152+
! input variables
153+
!-----------------------------------------------------------------
154+
155+
type (mpas_pool_type), intent(in) :: forcingPool
141156

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

159174
integer :: iEdge, kmin, kmax, k
160175
real (kind=RKIND) :: phase, u, v, time
176+
real (kind=RKIND), pointer :: forcingTimeIncrement
161177

162178
! End preamble
163179
!-----------------------------------------------------------------
164180
! Begin code
165181

166182
if (.not. config_use_manufactured_solution) return
167183

168-
time = daysSinceStartOfSim*86400.0_RKIND
184+
call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
185+
186+
time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement
169187

170188
do iEdge = 1, nEdgesOwned
171189

0 commit comments

Comments
 (0)