Merge branch 'jeremy-lilly/mpaso/fb-lts' into next (PR #6224)

Implement FB-LTS Time-Stepping Scheme

This PR implements a new local time-stepping (LTS) scheme in MPAS-Ocean
for single layer configurations. This new scheme, called FB-LTS, is
similar to to the LTS3 scheme already implemented by MPAS-O [1], but
boasts an increased maximal CFL number. This results in FB-LTS being
able to use time-steps between ~2 times larger than time-steps for LTS3
while costing the same number of tendency evaluations.

The FB-LTS scheme is based on the forward-backward, three-stage, second-
order Runge-Kutta (FB-RK(3,2)) scheme. The existing LTS3 scheme is based
on a three-stage, third-order, strong stability preserving Runge-Kutta
scheme (SSPRK3), so we expect FB-LTS to outperform LTS3 in admittable
time-step by factors of between ~1.6-2.2.

In preliminary performance experiments, FB-LTS provides a 8x speedup
over RK4 (and 2x over LTS3) in the hurricane Sandy test case, on a
variable resolution mesh going to resolutions as high as 2km in Delaware
Bay.

[BFB] - only for mpaso-standalone

components/mpas-ocean/bld/build-namelist

@@ -557,9 +557,9 @@ add_default($nl, 'config_Redi_constant_kappa');
 add_default($nl, 'config_Redi_maximum_slope');
 add_default($nl, 'config_Redi_use_slope_taper');
 add_default($nl, 'config_Redi_use_surface_taper');
+add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_use_quasi_monotone_limiter');
 add_default($nl, 'config_Redi_quasi_monotone_safety_factor');
-add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_Redi_horizontal_taper');
 add_default($nl, 'config_Redi_horizontal_ramp_min');

components/mpas-ocean/bld/build-namelist-section

@@ -96,9 +96,9 @@ add_default($nl, 'config_Redi_constant_kappa');
 add_default($nl, 'config_Redi_maximum_slope');
 add_default($nl, 'config_Redi_use_slope_taper');
 add_default($nl, 'config_Redi_use_surface_taper');
+add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_use_quasi_monotone_limiter');
 add_default($nl, 'config_Redi_quasi_monotone_safety_factor');
-add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_Redi_horizontal_taper');
 add_default($nl, 'config_Redi_horizontal_ramp_min');

components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml

@@ -146,15 +146,14 @@
 <config_Redi_maximum_slope>0.01</config_Redi_maximum_slope>
 <config_Redi_use_slope_taper>.true.</config_Redi_use_slope_taper>
 <config_Redi_use_surface_taper>.true.</config_Redi_use_surface_taper>
+<config_Redi_limit_term1>.true.</config_Redi_limit_term1>
 <config_Redi_use_quasi_monotone_limiter>.true.</config_Redi_use_quasi_monotone_limiter>
 <config_Redi_quasi_monotone_safety_factor>0.9</config_Redi_quasi_monotone_safety_factor>
-<config_Redi_limit_term1>.true.</config_Redi_limit_term1>
 <config_Redi_min_layers_diag_terms>0</config_Redi_min_layers_diag_terms>
 <config_Redi_min_layers_diag_terms ocn_grid="ARRM10to60E2r1">15</config_Redi_min_layers_diag_terms>
 <config_Redi_horizontal_taper>'ramp'</config_Redi_horizontal_taper>
 <config_Redi_horizontal_taper ocn_grid="SOwISC12to60E2r4">'RossbyRadius'</config_Redi_horizontal_taper>
 <config_Redi_horizontal_taper ocn_grid="ECwISC30to60E2r1">'RossbyRadius'</config_Redi_horizontal_taper>
-<!-- To do: ramp for WC but RossbyRadius for Cryo -->
 <config_Redi_horizontal_taper ocn_grid="IcoswISC30E3r5">'ramp'</config_Redi_horizontal_taper>
 <config_Redi_horizontal_ramp_min>20e3</config_Redi_horizontal_ramp_min>
 <config_Redi_horizontal_ramp_min ocn_grid="WCAtl12to45E2r4">30e3</config_Redi_horizontal_ramp_min>
@@ -180,7 +179,6 @@
 <config_GM_closure ocn_grid="WCAtl12to45E2r4">'constant'</config_GM_closure>
 <config_GM_closure ocn_grid="SOwISC12to60E2r4">'N2_dependent'</config_GM_closure>
 <config_GM_closure ocn_grid="ECwISC30to60E2r1">'N2_dependent'</config_GM_closure>
-<!-- To do: constant for WC but N2_dependent for Cryo -->
 <config_GM_closure ocn_grid="IcoswISC30E3r5">'constant'</config_GM_closure>
 <config_GM_constant_kappa>900.0</config_GM_constant_kappa>
 <config_GM_constant_kappa ocn_forcing="datm_forced_restoring" ocn_grid="oEC60to30v3wLI">600.0</config_GM_constant_kappa>
@@ -198,7 +196,6 @@
 <config_GM_spatially_variable_baroclinic_mode>3.0</config_GM_spatially_variable_baroclinic_mode>
 <config_GM_spatially_variable_baroclinic_mode ocn_grid="SOwISC12to60E2r4">1.0</config_GM_spatially_variable_baroclinic_mode>
 <config_GM_spatially_variable_baroclinic_mode ocn_grid="ECwISC30to60E2r1">1.0</config_GM_spatially_variable_baroclinic_mode>
-<!-- To do: 3.0 for WC but 1.0 for Cryo? -->
 <config_GM_spatially_variable_baroclinic_mode ocn_grid="IcoswISC30E3r5">3.0</config_GM_spatially_variable_baroclinic_mode>
 <config_GM_Visbeck_alpha>0.13</config_GM_Visbeck_alpha>
 <config_GM_Visbeck_max_depth>1000.0</config_GM_Visbeck_max_depth>
@@ -209,7 +206,6 @@
 <config_GM_horizontal_taper>'ramp'</config_GM_horizontal_taper>
 <config_GM_horizontal_taper ocn_grid="SOwISC12to60E2r4">'RossbyRadius'</config_GM_horizontal_taper>
 <config_GM_horizontal_taper ocn_grid="ECwISC30to60E2r1">'RossbyRadius'</config_GM_horizontal_taper>
-<!-- To do: ramp for WC but RossbyRadius for Cryo -->
 <config_GM_horizontal_taper ocn_grid="IcoswISC30E3r5">'ramp'</config_GM_horizontal_taper>
 <config_GM_horizontal_ramp_min>20e3</config_GM_horizontal_ramp_min>
 <config_GM_horizontal_ramp_min ocn_grid="WCAtl12to45E2r4">30e3</config_GM_horizontal_ramp_min>

components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml

@@ -202,9 +202,9 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_time_integrator" type="char*1024"
 	category="time_integration" group="time_integration">
-Time integration method.
+Time integration method. These options are only supported in standalone, not E3SM: 'LTS', 'FB_LTS'.
 
-Valid values: 'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit', 'LTS', 'split_explicit_ab2'
+Valid values: 'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit', 'split_explicit_ab2', 'LTS', 'FB_LTS'
 Default: Defined in namelist_defaults.xml
 </entry>
 
@@ -2038,7 +2038,7 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_n_ts_iter" type="integer"
 	category="split_timestep_share" group="split_timestep_share">
-number of large iterations over stages 1-3
+number of large iterations over stages 1-3; For the split_explicit_ab2 time integrator, this value only affects the first time step when it is not a restart run. For restart runs, this value has no effect on the split_explicit_ab2 time integrator.
 
 Valid values: any positive integer, but typically 1, 2, or 3
 Default: Defined in namelist_defaults.xml

components/mpas-ocean/src/Registry.xml

@@ -201,8 +201,8 @@
 					possible_values="Any time stamp in 'YYYY-MM-DD_hh:mm:ss' format. Items can be removed from the left if they are unused."
 		/>
 		<nml_option name="config_time_integrator" type="character" default_value="split_explicit_ab2"
-					description="Time integration method."
-					possible_values="'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit', 'LTS', 'split_explicit_ab2'"
+					description="Time integration method. These options are only supported in standalone, not E3SM: 'LTS', 'FB_LTS'."
+					possible_values="'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit', 'split_explicit_ab2', 'LTS', 'FB_LTS'"
 		/>
 		<nml_option name="config_number_of_time_levels" type="integer" default_value="2"
 					description="The number of time levels in the time-stepping scheme. This is used for array allocation."
@@ -1196,6 +1196,20 @@
 					possible_values="Any positive integer greater than or equal to one. A value of one employs the same dt in all regions."
 		/>
 	</nml_record>
+    <nml_record name="forward_backward" mode="forward">
+		<nml_option name="config_fb_weight_1" type="real" default_value="0.531"
+					description="The forward-backward weight for the first stage of FB-RK(3,2), used in FB_LTS."
+					possible_values="Any positive real number less than or equal to one."
+		/>
+		<nml_option name="config_fb_weight_2" type="real" default_value="0.531"
+					description="The forward-backward weight for the second stage of FB-RK(3,2), used in FB_LTS."
+					possible_values="Any positive real number less than or equal to one."
+		/>
+		<nml_option name="config_fb_weight_3" type="real" default_value="0.313"
+					description="The forward-backward weight for the third stage of FB-RK(3,2), used in FB_LTS."
+					possible_values="Any positive real number less than or equal to one."
+		/>
+	</nml_record>
 	<nml_record name="pressure_gradient" mode="forward">
 		<nml_option name="config_pressure_gradient_type" type="character" default_value="pressure_and_zmid"
 					description="Form of pressure gradient terms in momentum equation. For most applications, the gradient of pressure and layer mid-depth are appropriate.  For isopycnal coordinates, one may use the gradient of the Montgomery potential.  The sea surface height gradient (ssh_gradient) option is for barotropic, depth-averaged pressure."

components/mpas-ocean/src/mode_forward/Makefile

@@ -6,13 +6,19 @@ OBJS = mpas_ocn_forward_mode.o \
        mpas_ocn_time_integration_si.o \
        mpas_ocn_time_integration_split.o \
        mpas_ocn_time_integration_lts.o \
+       mpas_ocn_time_integration_fblts.o \
        mpas_ocn_time_integration_split_ab2.o
 
 all: forward_mode
 
 forward_mode: $(OBJS) 
 
-mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o mpas_ocn_time_integration_si.o  mpas_ocn_time_integration_split.o mpas_ocn_time_integration_lts.o mpas_ocn_time_integration_split_ab2.o
+mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o \
+	                         mpas_ocn_time_integration_si.o \
+							 mpas_ocn_time_integration_split.o \
+							 mpas_ocn_time_integration_lts.o \
+							 mpas_ocn_time_integration_fblts.o \
+							 mpas_ocn_time_integration_split_ab2.o
 
 mpas_ocn_time_integration_rk4.o: 
 
@@ -27,6 +33,7 @@ mpas_ocn_forward_mode.o: mpas_ocn_time_integration.o \
                          mpas_ocn_time_integration_si.o \
                          mpas_ocn_time_integration_split.o \
                          mpas_ocn_time_integration_lts.o \
+					     mpas_ocn_time_integration_fblts.o \
                          mpas_ocn_time_integration_split_ab2.o
 
 clean:

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

@@ -34,6 +34,7 @@ module ocn_time_integration
    use ocn_time_integration_split
    use ocn_time_integration_si
    use ocn_time_integration_lts
+   use ocn_time_integration_fblts
    use ocn_time_integration_split_ab2
 
 
@@ -72,7 +73,8 @@ module ocn_time_integration
        timeIntSemiImplicit     = 3, &! Semi-implicit
        timeIntRK4              = 4, &! 4th-order Runge-Kutta
        timeIntLTS              = 5, &! local time-stepping
-       timeIntSplitExplicitAB2 = 6   ! split-explicit AB2 baroclinic
+       timeIntFBLTS            = 6, &! forward-backward lts
+       timeIntSplitExplicitAB2 = 7   ! split-explicit AB2 baroclinic
 
 !***********************************************************************
 
@@ -136,6 +138,8 @@ subroutine ocn_timestep(domain, dt, timeStamp)!{{{
          call ocn_time_integrator_rk4(domain, dt)
       case (timeIntLTS)
          call ocn_time_integrator_lts(domain, dt)
+      case (timeIntFBLTS)
+         call ocn_time_integrator_fblts(domain, dt)
       case (timeIntSplitExplicitAB2)
          call ocn_time_integrator_split_ab2(domain, dt)
       end select
@@ -232,6 +236,10 @@ subroutine ocn_timestep_init(domain, dt, err)!{{{
       case ('LTS')
          timeIntegratorChoice = timeIntLTS
          call ocn_time_integration_lts_init(domain)
+
+      case ('FB_LTS')
+         timeIntegratorChoice = timeIntFBLTS
+         call ocn_time_integration_fblts_init(domain)
 
       case ('split_explicit_ab2')
          timeIntegratorChoice = timeIntSplitExplicitAB2